/* *----------------------------------------------------------------------* * * * CERN * * * * European Organization for Nuclear Research * * * * GSI * * * * Gesellschaft fuer Schwerionenforschung * * * * Program name: HEADTAIL Version 2.0 * * * * Author and contact: G. RUMOLO * * Hochstrom-Strahlphysik * * BR1 3.011 * * GSI, Planckstrasse 1 * * D-64291 Darmstadt * * GERMANY * * * * contact: F. ZIMMERMANN * * AB (formerly) SL Division * * CERN * * CH-1211 GENEVA 23 * * SWITZERLAND * * * * Tel. [041] (022) 767 9054 * * frank.zimmermann@cern.ch * * [049] 6159 71-298523 * * grumolo@gsi.de * * * * Copyright CERN, Geneva 2003 - Copyright and any other * * appropriate legal protection of this computer program and * * associated documentation reserved in all countries of the * * world. * * * * Organizations collaborating with CERN may receive this program * * and documentation freely and without charge. * * * * CERN undertakes no obligation for the maintenance of this * * program, nor responsibility for its correctness, and accepts * * no liability whatsoever resulting from its use. * * * * Program and documentation are provided solely for the use of * * the organization to which they are distributed. * * * * This program may not be copied or otherwise distributed * * without permission. This message must be retained on this and * * any other authorized copies. * * * * The material cannot be sold. CERN should be given credit in * * all references. * * * *----------------------------------------------------------------------* In this version the possibility of off-centering the beam in t=0 has also been introduced. There are two more parameters to be input, xkick and ykick, which define the kick amplitude. The tune spread introduced by space charge has also been updated: the tune shift of each particle is calculated at each turn taking into account the local and instantaneous sigma in the tune shift Laslett formula. Actually, according to the switch i_rot, it's possible to get the space charge roatation at each turn either around the local centroids and by using the local sigmas (i_rot=0), or around the full bunch centroid and by using the total sigmas (i_rot=1). i_rot is the second last entry of the new input file. The option of the electrons moving in a sinusoidal field is here included - it corresponds to i_dip=2 and takes the input value for the solenoid field (bsol, last line in the input file) The option of the proton space charge has been added - by using the routine ffrank, we make each bunch slice act on each proton in the slice, too. A new switch has been introduced, i_space, which is a value to be input in the input file. The most updated input file is to be found in the directory ~/Elcloud/HEADTAIL/DATA/PIC/TUNESPR and it's nobb0quacin.cfg - sisexp.cfg is in the ../PIC/EXPER directory (input values: 37) !!!! This file has two executables: !!!! xcloud: electron-cloud problems oriented - it has NPR=300000 NEL=100000 NBIN=50 ximped: impedance problems oriented - it has NPR=1000000 NEL=1 NBIN=70 Routine new_fraction is introduced, in which 500 macroelectrons are newly initialized at each interaction slice-cloud. This is done only when a solenoid field is acting. It models the fact that electrons in the solenoid field are actually moving longitudinally. We take these numbers: electron cooler section 3m long, bunch 5m long, electron velocity is 0.145c and proton (or D+) velocity 0.947c. We easily find Deltat_b = 5/0.947c = 17.6ns -> during one full bunch pass the electrons will have advanced by Deltas_e = 0.145c Deltat_b = 0.765m. Therefore 1/4 of the total cooler length. This means that at the end of the interaction, 1/4 of the cloud keeps actually a pregressively "shorter memory" of the previous interaction. If after 50 slices 1/4 of the macroelectrons (25000) will have this partial memory, the effect can be achieved if at each slice we regenerate 25000/50 = 500 electrons. Last improvement has been made with Daniel Schulte, he has basically made the PIC field calculation faster by changing the modules fourtrans with improved versions (fourtrans2 and fourtrans3) which link to the file fourtrans2.c, included in the present program, and to routines from the fftw library which has also been copied in the fftw directory from my home directory (7/12/2001). As the amplitude detuning has been also introduced, there is a new switch i_oct which sets it working or not. It is the last entry in the input file (line 38). If amplitude detuning is turned on, its parametrization is given in the 'init_values' routine of this code and actually reflects the parametrization expected for SPS (from Rogelio's note). 11/04/2002 We have introduced a new option: i_dip=3. The electrons move in the strong combined functions magnetic field, like in PS. With the assumption of strong field, we have simply confined the electron motion on trajectories that follow the field lines. 29/04/2002 The new flag i_incoh has been introduced to study the instability free of the monopole mode as suggested by Evgueni Perevedentsev. The bunch phase space centroid is subtracted turn by turn and slice by slice from the particle phase space coordinates to concentrate only on the e-cloud induced emittance growth The flag elflag defines the electron distribution: 1- Rectangular 2- Elliptical 3- Stripe-like rectangular 4- Parabolic Number of entries of the input file: 40 !!! new entries to be added: Coherent_centroid_motion_(0->off_1->on): Electron_distribution(1->Rect_2->Ellipt_3->Stripe_4->Parab): for most updated input file look into: /afs/cern.ch/user/g/grumolo/Elcloud/HEADTAIL/DATA/PIC/PEREVED/eugeni*.cfg /afs/cern.ch/user/g/grumolo/Elcloud/HEADTAIL/DATA/PIC/PS/COMBIN/cfcurr8ott.cfg 3/06/2002 new flag i_coupl to have possibly linear xy-coupling in the machine. The parameter key_xy is also defined to tune the strength of the coupling. To be input on last line of the input file. Number of entries of the input file: 42 !!! new entries to be added: Linear_coupling_switch(1->on_0->off): 1 Linear_coupling_coefficient_[1/m]: 0.015 17/06/2002 ->new<- to take into account horizontal dispersion: the value is input in the variable dispx (43rd line in the input file !!), and the contribution of dispersion comes in when the interaction between bunch and cloud is calculated. This introduces a x-z correlation first, and obviously it changes the oscillation frequency of the electrons because it changes the horizontal rms-size value. Also the 2-stripes distribution is included for electrons, which takes into account how electrons are actually distributed inside a dipole field. Number of entries of the input file: 45 !!! new entries: Average_dispersion_function_in_the_ring_[m]: 2.0 Position_of_the_stripes_[units_of_sigmax]: 3.0 Width_of_the_stripes_[units_of_sigmax]: 0.5 and the line on the electron distribution becomes: Electron_distrib_(1->Rect_2->Ellip_3->[1_strp]_4->[2_strp]_5->Parab): 08/08/2002 A kick in the longitudinal direction has been introduced, after having changed the bining in such a way that the code can deal with bunch dipole oscillations in z and variations of sigma_z during the evolution (to model for example a debunching beam). New line in the input file is therefore needed: Kick_in_the_longitudinal_direction_[m]: 0.03 If isyn==2 we also give the possibility to the bunch to de-bunch. Synchrotron motion is off and the particles follow their longitudinal straight trajectories according to their initial momenta. 13/08/02 Changing the transverse impedance: the transverse force due to the broad band impedance has been adapted to a case of flat chamber (the two limit cases are now possible: round chamber i_pipe=0 (a=b) or very flat chamber i_pipe=1 (a>>b)) this introduces one more line in the input file to set i_pipe, but also two lines of the input need to be taken out because we have one single set of parameters for the trasverse impedance, which will be used both for horizontal and vertical kicks + the flat chamber correction. OLD: Switch_for_wake_fields: 0 Res_frequency_of_broad_band_resonator_[GHz]: 1.3 Horizontal_quality_factor: 1. Vertical_quality_factor: 1. Shunt_impedance_x_[MOhm/m]: 5. Shunt_impedance_y_[MOhm/m]: 10. NEW: Switch_for_wake_fields: 0 Switch_for_pipe_geometry_(0->round_1->flat): 1 Res_frequency_of_broad_band_resonator_[GHz]: 1.3 Transverse_quality_factor: 1. Transverse_shunt_impedance_[MOhm/m]: 10. 15/08/02 It is possible to store the bunch shape every n_diag turns. This makes sense when we are in debunching mode or in bunched mode with a bunch not well matched to the bucket. One more input parameter: Number_of_turns_between_two_bunch_shape_acquisitions: 20 Longitudinal wake field has been introduced due to a resonator. We want to model the parasitic effect of a cavity on a single bunch. We can neglect the turn-to-turn effect, since the quality factor is such (~150) that the field has completely decayed when the bunch goes again through the same cavity where it was generated (Q=5000 is at the limit if this approximation for a frequency around 200 MHz). Again we need to input three numbers to model this resonator: Res_frequency_of_broad_band_resonator_[GHz]: 1.3 Transverse_quality_factor: 1. Transverse_shunt_impedance_[MOhm/m]: 10. Res_frequency_of_longitudinal_resonator_[MHz]: 200 Longitudinal_quality_factor: 140. Longitudinal_shunt_impedance_[MOhm]: 5. 20/08/02 Fields with conducting boundary conditions on a rectangular box can be calculated instead of the simple periodic condition in x and y. The method consists in unfolding the charge distribution on a double-sized grid on which the image charges are placed. The routine image_unfold has been introduced to accomplish this. One more input needed: Boundary_condition_(0->Periodic_1->Conducting): 1 */ /* Include all facilities */ #include #include //#include #include #include #include #include #include /* with this package one can define complex variables, if needed; for instance, complex z */ /* Definition of constants */ /* Physical constants.... */ #define PI 3.14159265358979 #define C 2.99792458e8 /* velocity of light [m/s] */ #define E 1.6021892e-19 /* electron charge [Cb] */ #define MU 1.66055e-27 /* proton mass [kg] */ #define ME 9.109534e-31 /* proton-positron mass [kg] */ #define EPS0 8.8542e-12 /* dielectric constant [C^2/Jm] */ #define RP 1.5458775e-18 /* classical proton radius [m] */ #define RE 2.817938e-15 /* classical electron radius [m] */ #define I0 31.356e6 /* ... for the tune shift.. [A] */ #define ELMS 9.109534e-31 /* electron mass [Kg] */ /* Parameters for the computation... */ #define NPR 1000000 /* Number of macroparticles (bunch) to be used */ #define NPR1 90000 /* Maximum number of macroparticles in a slice */ #define NEL 1 /* Number of macroelectrons (cloud) to be used */ #define NBIN 200 /* Number of 'slices' in which each single bunch is sub-divided */ #define NDC 128 /* Number of cells to produce a distribution function */ /* type definition */ typedef struct /* GRID */ { int n_cell_x,n_cell_y; int timestep; double delta_x,delta_y; double offset_x,offset_y; double min_x,max_x,min_y,max_y; double cut_x,cut_y,cut_z; double step,n_macro_1,n_macro_2,step1,step2,scal_step[2]; float rho_x_1,rho_y_1,rho_sum_1,rho_x_2,rho_y_2,rho_sum_2; double rho_factor; int n_m_1,n_m_2; double *rho1,*rho2,*phi1,*phi2,*dist,*temp; int integration_method,distribution,interpolation; } GRID; /* Definition of the external variables */ FILE *bunch_pr , *headtail_pr , *prot_phase , *prot_inphase , *samp_prot , *prot_distr ; char filename[80] ; long iind[NBIN] , /* iind[i] gives the number of macroparticles in the i-th slice of bunch */ indexp[NBIN][NPR1] , /* in the bunch slicing, indexp[i][j] contains the original bunch index of the particle that has index j in the i-th slice */ imain , /* down to mmain -> indices used in main for loops */ jmain , kmain , lmain , mmain , elflag , /* flag for the e-distribution: see above */ it , /* index used for the main cycle over the total number of steps: turns * kicks */ ind , nturn , /* number of simulated turns -> input parameter */ nkick , /* number of kicks per turn -> input parameter */ nstep , /* total number of steps: nturn * nkick */ isyn , /* flag for the synchrotron motion: input parameter */ ipr_pos , /* flag for proton/positron bunch: input parameter */ imesize , /* scale factor for macroelectron size: input parameter */ ielsizex , /* multiplication factor for the e- cloud x-size: input parameter */ ielsizey , /* multiplication factor for the e- cloud y-size: input parameter */ ili , /* flag for the function ffrank: it has to be =! 1 */ iwake_flag , /* flag for the wake field: input parameter */ ishift , /* flag to activate the tune shift option */ sl_ind , /* it is used as slice index when the incoherent space charge tune shift is added */ i_method , /* flag to choose the field calculation method */ i_rot , /* switch for the space charge tune rotation */ i_space , /* Flag for the space charge forces */ i_oct , /* Flag for the action dependent tunes */ i_dip , /* flag to switch off the electron horizontal motion in order to simulate their motion in a strong vertical dipole field */ i_incoh , /* flag to investigate on the monopole mode as suggested by Evgueni Perevedentsev */ i_coupl , /* flag to switch on and off linear coupling */ i_pipe , /* flag for pipe shape in effective wake calculation 0 -> circular (a=b) 1 -> flat (a>>b) */ n_diag , /* every how many turns the bunch shape is stored */ hh , /* Harmonic number of sinusoidal bucket */ npr_eff ; /* number of bunch particles lost turn by turn */ double xpr[NPR] , /* horizontal coordinate for protons or positrons */ ypr[NPR] , /* vertical coordinate for protons or positrons */ xpp[NPR] , /* dx/ds for protons or positrons */ ypp[NPR] , /* dy/ds for protons or positrons */ sz[NPR] , /* longitudinal coordinate for protons or positrons */ dp[NPR] , /* relative momentum spread for protons or positrons */ xe[NEL] , /* horizontal coordinate for electrons */ ye[NEL] , /* vertical coordinate for electrons */ xpe[NEL] , /* dxe/dt for electrons */ ype[NEL] , /* dye/dt for electrons */ //xe0[NEL] , //ye0[NEL] , zs[NBIN] , /* longitudinal position of each bunch slice */ npr[NBIN] , /* number of particles in each slice */ sxs[NBIN] , /* horizontal rms value of each slice */ sys[NBIN] , /* vertical rms value of each slice */ xs[NBIN] , /* horizontal mean offset of each slice */ ys[NBIN] , /* vertical mean offset of each slice */ xps[NBIN] , /* horizontal mean dx/ds of each slice */ yps[NBIN] , /* vertical mean dy/ds of each slice */ indx[NDC] , /* grid points for the horizontal distribution of electrons */ //hde[NDC] , /* horizontal distribution of electrons after having been kicked */ //indy[NDC] , /* grid points for the vertical distribution of electrons */ //vde[NDC] , /* vertical distribution of electrons after having been kicked */ indxp[NDC] , /* grid points for the bunch long. distribution */ hdp[NDC] , /* Bunch longitudinal distribution */ //indyp[NDC] , /* grid points for the vertical distribution of the bunch particles */ //vdp[NDC] , /* vertical distribution of bunch particles before interacting with the particles */ //indxpp[NDC] , /* grid points for the distribution of horizontal momenta of the bunch particles */ //hpdp[NDC] , /* distribution of horizontal momenta of the bunch particles before interacting with the electrons */ //indypp[NDC] , /* grid points for the distribution of vertical momenta of the bunch particles */ //vpdp[NDC] , /* distribution of vertical momenta of the bunch particles before interacting with the electrons */ vprx[NPR1] , /* temporary vector that contains all the x's of the bunch particles in a slice */ vpry[NPR1] , /* temporary vector that contains all the y's of the bunch particles in a slice */ vprxp[NPR1] , /* temporary vector that contains all the xp's of the bunch particles in a slice */ vpryp[NPR1] , /* temporary vector that contains all the yp's of the bunch particles in a slice */ /* now the declaration of the input parameters starts.... it goes up to nus = synchrotron tune */ sz0 , /* [m] initial rms-bunch length */ dp0 , /* initial rms-momentum spread */ sx0 , /* [m] initial rms hor. beam size */ sx0_cd , /* [m] rms hor. beam size dispersion corrected */ sy0 , /* [m] initial rms vert. beam size */ npr0 , /* number of particles in the bunch */ betax0 , /* [m] average hor. beta function around the ring */ betay0 , /* [m] average vert. beta function around the ring */ dispx , /* [m] average dispersion function in the ring */ dcloud , /* average number of electrons per m^3 in the ring */ circ , /* [m] ring circumference */ gammar , /* relativistic gamma of the bunch particles */ htune , /* horizontal tune of the machine */ vtune , /* vertical tune of the machine */ csix , /* horizontal chromaticity of the machine */ csiy , /* vertical chromaticity of the machine */ alpha , /* momentum compaction factor of the machine */ nus , /* synchrotron tune */ xkick , /* [m] Amplitude of the horizontal kick in t=0 */ ykick , /* [m] Amplitude of the vertical kick in t=0 */ bsol , /* [T] Strength of the solenoid field */ bdip , /* [T] Strength of the dipole field in a combined functions magnet (like in the PS) */ bgrad , /* [T/m] Gradient of the magnetic field in a combined functions magnet (like in the PS) */ rpipex , /* horizontal semiaxis of the beam pipe: gives a x-size to the e- cloud */ rpipey , /* vertical semiaxis of the beam pipe: gives a y-size to the e- cloud */ gridextx , /* Horizontal extension of the grid for the fields PIC calculation */ gridexty , /* Vertical extension of the grid for the fields PIC calculation */ strpos , /* position of the stripes in units of sigma_x */ strsig , /* width of the stripes in units of sigma_x */ eta , /* slip factor: alpha - 1/gammar^2 */ omega0 , /* bunch revolution frequency */ omegas , /* synchrotron frequency */ factor1 , /* auxiliary parameter for the force calculation */ factor2 , /* auxiliary parameter for the force calculation */ fscale1 , /* auxiliary parameter for the force calculation */ fscale2 , /* auxiliary parameter for the force calculation */ eel , /* auxiliary parameter for the force calculation */ epr , /* auxiliary parameter for the force calculation */ wakefac , /* factor for the wake field kick */ tstep , /* time step during the interaction bunch-cloud */ dlstep , /* s-distance between two kick points */ tstepp , /* time step between two kick points */ qp0 , /* number of bunch particles in one macroparticle */ qe0 , /* number of electrons in one macroelectron */ p0 , /* average momentum of the bunch particles */ sx2 , /* finite macroelectron horizontal size */ sy2 , /* finite macroelectron vertical size */ key_xy , /* [1/m] strength of linear coupling */ cxa , /* cxa -> ssa: transport parameters for the phase */ sxa , /* space coordinates of the bunch particles (linear */ cya , /* transport + chromaticity kick in the transverse */ sya , /* plane, linear synchrotron motion longitudinally) */ csa , ssa , dtunxmax , /* maximum space charge horizontal tune shift in the bunch */ dtunymax , /* maximum space charge vertical tune shift in the bunch */ dtunx , dtuny , lambda , /* local variable to store the exponential factor for the tune shift modulation along the bunch */ cent_x , /* cent_*: local variables to store the centroid */ cent_y , /* position and slope of each slice used for the */ cent_xp , /* incoherent space charge tune spread */ cent_yp , xx1 , /* xx1->yyp1: local variables used in the particles */ xxp1 , /* transport around the ring */ yy1 , yyp1 , epsx , /* single particle horizontal emittance variable */ epsy , /* single particle vertical emittance variable */ qpxx , /* Q_x = Q_x0 + eps_x * qpxx + eps_y * qpxy */ qpyy , /* Q_y = Q_y0 + eps_x * qpxy + eps_y * qpyy */ qpxy , bin_frac , /* fraction of bunch particles that are binned */ xave , /* average horizontal offset of the beam */ yave , /* average vertical offset of the beam */ zave , /* average longitudinal offset of the beam */ dpmean , /* average long. momentum offset of the beam */ xpave , /* average beam offset in dx/ds */ ypave , /* average beam offset in dy/ds */ sxsize , /* square of horizontal rms bunch size */ sysize , /* square of vertical rms bunch size */ szsize , /* longitudinal rms bunch size */ sxpsize , sypsize , timep , /* time at which the kick occurs */ xemit , /* horizontal emittance */ yemit , /* vertical emittance */ invarx , /* Courant-Snyder horizontal invariant */ invary , /* Courant-Snyder vertical invariant */ corryz , /* y-z correlation */ npr1 , /* number of protons contained in one slice */ xoff , /* slice horizontal off-center */ yoff , /* slice vertical off-center */ sx , /* slice horizontal rms size */ sy , /* slice vertical rms-size */ xta , yta , zstep , /* length of a bunch slice */ zi , /* longitudinal negative coordinate (head of the bunch is zi = 0) for wake field */ x1 , /* macroelectron x position at which the force from each slice is calculated */ ips1 , /* macroelectron y position at which the force from each slice is calculated */ pre , /* coefficient for the y-component of the force */ pim , /* coefficient for the x-component of the force */ kick_xe , /* velocity x-kick from each slice on each macroelectron */ kick_ye , /* velocity y-kick from each slice on each macroelectron */ xof2 , /* macroelectron x-position */ yof2 , /* macroelectron y-position */ x2 , /* x-position of the bunch particle on which the macroelectron acts */ y2 , /* y-position of the bunch particle on which the macroelectron acts */ xof3 , /* initial x-position of the macroelectron - used to clean out fluctuations */ yof3 , /* initial y-position of the macroelectron - used to clean out fluctuations */ xta2 , yta2 , pre2 , /* coefficient for the x-component of the force */ pim2 , /* coefficient for the y-component of the force */ kick_x , /* x-kick on each bunch particle */ kick_y , /* y-kick on each bunch particle */ kick_z , /* z-kick from on each bunch particle */ z0_kick , /* z-kick from each macroelectron on each bunch particle */ kick_xp_sc , /* dx/ds-kick from each bunch slice on each bunch particle */ kick_yp_sc , /* dy/ds-kick from each bunch slice on each bunch particle */ x0 , /* x-coordinate of a bunch particle after having been kicked, before transport to the next kick */ xp0 , /* dx/ds of a bunch particle after having been kicked and before transport to the next kick */ ips0 , /* y-coordinate of a bunch particle after having been kicked, before transport to the next kick */ yp0 , /* dy/ds of a bunch particle after having been kicked and before transport */ szz0 , /* longitudinal coordinate of a bunch particle after having been kicked, before transport to the next kick */ dpp0 , /* momentum spread of a bunch particle after having been kicked, before transport to the next kick */ r11 , /* r11 -> r44, coefficients of the transport matrix for the transverse phase space coordinates between two kicks */ r12 , r21 , r22 , r33 , r34 , r43 , r44 , merit , /* quality factor of the broad band impedance */ freqr , /* resonant frequency of the broad band resonator to be expressed in GHz (x and y assumed equal) */ zetat , /* shunt impedance of the braod band resonator to be expressed in MOhm/m */ meritz , /* quality factor of the longitudinal impedance */ freqrz , /* resonant frequency of the longitudinal resonator to be expressed in MHz */ rshz , /* shunt impedance of the braod band resonator to be expressed in MOhm */ omegar , /* auxiliary parameters for transverse wake */ omegabar , alpha , omegarz , /* auxiliary parameters for longitudinal wake */ omegabarz , alphaz , voltm , /* Maximum voltage in sinusoidal bucket */ wrf , /* Radio-frequency of sinusoidal bucket */ coeff3 , /* Auxiliary parameters for longitudinal motion */ coeff4 , /* in sinusoidal bucket */ coeff5 , coeffspch , /* coefficient for longitudinal space charge in a sinusoidal bucket */ k1 , k2 , k3 , k4 , kp1 , kp2 , kp3 , kp4 ; /* For sake of clarity we specify here that the broad band model that we've assumed considers an impedance (omegar = 2*PI*freqr*10^(-9)): Z = omegar/omega * zetat/[1+i merit (omegar/omega - omega/omegar)]. If the impedance is: Z = c/omega * R_s/[1+i merit (omegar/omega - omega/omegar)], then one has to convert R_s into our input parameter by means of: zetat=c/omegar*R_s zetat [Ohm/m] and R_s [Ohm/m^2] In the longitudinal plane: Z|| = rshz/[1+i meritz (omegarz/omega - omega/omegarz)] */ /* Global declarations... function that evaluates the electric mutual force between bunch particles and electrons has to be included */ extern void ffrank (double*, double*, double*, double*, double*, double*, double*, double*, long*); /****************************************************************************** *** Subroutine : read_data *** *** Effect : Define main parameters by reading data from specified *** *** configuration-file *** *** Parameters : none *** *** Gbl var used : *** *** Gbl var effect : *** *** Constants used : none *** *** Subrout. used : none *** ******************************************************************************/ read_data () { char cfg_filename[30], dummy_string[64]; FILE *cfg_file_ptr; printf ("\n\n Please specify name (without extention) of desired") ; printf ("\n configuration-file : ") ; scanf ("%s",filename) ; strcpy (cfg_filename,filename); strcat (cfg_filename,".cfg"); cfg_file_ptr = fopen (cfg_filename,"r"); if (cfg_file_ptr == NULL) { printf ("\n\n Configuration-file does not exist.\n") ; end () ; } fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&ipr_pos); //printf ("ipr_pos=%d\n",ipr_pos); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&dcloud); //printf ("dcloud=%lf\n",dcloud); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&npr0); //printf ("npr0=%lf\n",npr0); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&betax0); //printf ("betax0=%lf\n",betax0); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&betay0); //printf ("betay0=%lf\n",betay0); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&sz0); //printf ("sz0=%lf\n",sz0); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&sx0); //printf ("sx0=%lf\n",sx0); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&sy0); //printf ("sy0=%lf\n",sy0); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&dp0); //printf ("dp0=%lf\n",dp0); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&nus); //printf ("nus=%lf\n",nus); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&alpha); //printf ("alpha=%lf\n",alpha); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&circ); //printf ("circ=%lf\n",circ); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&gammar); //printf ("gammar=%lf\n",gammar); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&nkick); //printf ("nkick=%d\n",nkick); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&nturn); //printf ("nturn=%d\n",nturn); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&ielsizex); //printf ("ielsizex=%d\n",ielsizex); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&ielsizey); //printf ("ielsizey=%d\n",ielsizey); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&htune); //printf ("htune=%lf\n",htune); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&vtune); //printf ("vtune=%lf\n",vtune); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&csix); //printf ("csix=%lf\n",csix); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&csiy); //printf ("csiy=%lf\n",csiy); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&isyn); //printf ("isyn=%d\n",isyn); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&imesize); //printf ("imesize=%d\n",imesize); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&iwake_flag); //printf ("iwake_flag=%d\n",iwake_flag); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&i_pipe); //printf ("i_pipe=%d\n",i_pipe); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&freqr); //printf ("freqr=%lf\n",freqr); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&merit); //printf ("merit=%lf\n",merit); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&zetat); //printf ("zetat=%lf\n",zetat); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&freqrz); //printf ("freqrz=%lf\n",freqrz); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&meritz); //printf ("meritz=%lf\n",meritz); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&rshz); //printf ("rshz=%lf\n",rshz); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&ishift); //printf ("ishift=%d\n",ishift); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&i_method); //printf ("i_method=%d\n",i_method); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&i_dip); //printf ("i_dip=%d\n",i_dip); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&xkick); //printf ("xkick=%lf\n",xkick); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&ykick); //printf ("ykick=%lf\n",ykick); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&i_space); //printf ("i_space=%d\n",i_space); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&i_rot); //printf ("i_rot=%d\n",i_rot); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&bsol); //printf ("bsol=%lf\n",bsol); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&i_oct); //printf ("i_oct=%d\n",i_oct); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&i_incoh); //printf ("i_incoh=%d\n",i_incoh); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&elflag); //printf ("elflag=%d\n",elflag); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&i_coupl); //printf ("i_coupl=%d\n",i_coupl); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&key_xy); //printf ("key_xy=%lf\n",key_xy); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&dispx); //printf ("dispx=%lf\n",dispx); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&strpos); //printf ("strpos=%lf\n",strpos); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&strsig); //printf ("strsig=%lf\n",strsig); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&z0_kick); //printf ("z0_kick=%lf\n",z0_kick); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&n_diag); //printf ("n_diag=%d\n",n_diag); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%lf",&voltm); //printf ("z0_kick=%lf\n",voltm); fscanf (cfg_file_ptr,"%s",dummy_string); fscanf (cfg_file_ptr,"%d",&hh); //printf ("n_diag=%d\n",hh); if(i_dip==2) { elflag = 2; printf("elflag set to 2, E-cooler case considered"); } printf ("\n\n\t\tConfiguration-file loaded successfully.\n") ; } /****************************************************************************** *** Subroutine : init_values *** *** Effect : Initialize all parameter values, necessary for *** *** calculation. *** *** Parameters : none *** *** Gbl var used : *** *** Gbl var effect : *** *** Constants used : *** *** Subrout. used : none *** ******************************************************************************/ init_values () { double area; eta = alpha - 1/gammar/gammar; omega0 = C/circ*2*PI; omegas = nus*omega0; tstep = 4.*sz0/(double)NBIN/C; factor1 = RE*C*sqrt(2*PI)/tstep; fscale1 = factor1/sqrt(2*PI); dlstep = circ/(double)nkick; tstepp = dlstep/C; npr_eff = 0; omegar = 2*PI*freqr*1e9; alpha = omegar/(2*merit); omegabar = sqrt(omegar*omegar - alpha*alpha); omegarz = 2*PI*freqrz*1e6; alphaz = omegarz/(2*meritz); omegabarz = sqrt(omegarz*omegarz - alphaz*alphaz); wrf = (double)hh*omega0; switch(ipr_pos){ case 1: factor2 = RP*sqrt(2*PI)/gammar/tstep; p0 = MU*gammar*C; break; case 2: factor2 = RE*sqrt(2*PI)/gammar/tstep; p0 = ME*gammar*C; break; case 3: factor2 = RP/2*sqrt(2*PI)/gammar/tstep; p0 = 2*MU*gammar*C; break; } coeff3 = eta/fabs(eta)*E*voltm/(p0*circ); coeff4 = wrf/C; coeff5 = -eta*C; coeffspch = (0.67 + 2*log(((double)ielsizex+(double)ielsizey)/2))*E*E*npr0/(MU*sqrt(2*PI)*4*PI*EPS0*gammar*gammar*p0); fscale2 = factor2/sqrt(2*PI); sx0_cd = sqrt(sx0*sx0 + dispx*dispx*dp0*dp0); rpipex = sx0_cd*(double)ielsizex; rpipey = sy0*(double)ielsizey; qp0 = npr0/(double)NPR; switch(elflag){ case 1: qe0 = 4.0*rpipex*rpipey*dlstep*dcloud/(double)NEL; break; case 2: qe0 = PI*rpipex*rpipey*dlstep*dcloud/(double)NEL; break; case 3: qe0 = 8.0*sx0_cd*rpipey*dlstep*dcloud/(double)NEL; break; case 4: qe0 = 2.0*rpipey*8.0*strsig*sx0*dlstep*dcloud/(double)NEL; break; case 5: qe0 = 4.0*rpipex*rpipey*dlstep*dcloud/(double)NEL; break; } area = qe0/dlstep/dcloud; sx2 = 0.5*sqrt(area*sx0_cd/sy0)/(double)imesize; sy2 = sy0/sx0_cd*sx2; cxa = cos(2*PI*htune/(double)nkick); sxa = sin(2*PI*htune/(double)nkick); cya = cos(2*PI*vtune/(double)nkick); sya = sin(2*PI*vtune/(double)nkick); dtunxmax = C*npr0*E*circ*circ/(4*PI*PI*I0*gammar*gammar*gammar*htune*sx0_cd*sx0_cd*sqrt(2*PI)*sz0); dtunymax = C*npr0*E*circ*circ/(4*PI*PI*I0*gammar*gammar*gammar*vtune*sy0*sy0*sqrt(2*PI)*sz0); csa = cos(omegas*tstepp); ssa = sin(omegas*tstepp); nstep = nkick*nturn; switch(ipr_pos){ case 1: wakefac = - E*E/(MU*gammar*C*C); break; case 2: wakefac = - E*E/(ME*gammar*C*C); break; case 3: wakefac = - E*E/(2*MU*gammar*C*C); break; } eel = 1.0; epr = 1.0; //i_oct = 1; /* SPS octupole setting: radial -0.875 and vertical -0.8 - signs have been reversed to check that now space charge should help the decoherence qpxx = 424.9; qpxy = -878; qpyy = 1155; */ /* Rogelio: SPS octupole setting: radial -0.875 and vertical -0.8 */ qpxx = -424.9; qpxy = 878; qpyy = -1155; /* Rogelio: SPS octupole setting: radial -0.875 and vertical 1.0 qpxx = -283.25; qpxy = 113.3; qpyy = -420; */ /* Coefficients of the magnetic field in the PS: B_x = - bgrad*y, B_y = bdip + bgrad*x The signs of the gradients do not affect the induced wake */ bdip = 1.2; bgrad = 5.2; } /****************************************************************************** *** Subroutine : build_distrib *** *** Effect : Calculate arbitrary distribution of bunch. *** *** Parameters : vector containing the quantity whose distribution *** *** we want to plot, number of elements of this vector *** *** Gbl var effect : *** *** Constants used : *** *** Subrout. used : none *** ******************************************************************************/ void verteil(double data[], unsigned long nele, double distri[], unsigned long ncell, double index[]) { long j , channel1 , channel2 , numb ; double scaled , mean_data , mean_sqr_data , fraction1 , fraction2 , data_center, data_width ; mean_data = 0 ; mean_sqr_data = 0 ; for (j=0 ; j 0.5) { fraction1 = 1.5 - scaled + channel1 ; channel2 = channel1 + 1 ; fraction2 = 1.0 - fraction1 ; } else { fraction1 = scaled - channel1 + 0.5 ; channel2 = channel1 - 1 ; fraction2 = 1.0 - fraction1 ; } if (channel1 >= 0 && channel1 < ncell) distri[channel1] += fraction1 ; if (channel2 >= 0 && channel2 < ncell) distri[channel2] += fraction2 ; } for (j=0 ; j 0.5) { fraction1 = 1.5 - scaled + channel1 ; channel2 = channel1 + 1 ; fraction2 = 1.0 - fraction1 ; } else { fraction1 = scaled - channel1 + 0.5 ; channel2 = channel1 - 1 ; fraction2 = 1.0 - fraction1 ; } if (channel1 >= 0 && channel1 < ncell) distri[channel1] += fraction1 ; if (channel2 >= 0 && channel2 < ncell) distri[channel2] += fraction2 ; } for (j=0 ; j= 1); xpe[i] = 0.0; ype[i] = 0.0; } if(elflag==3) for (i=0; i 0.5) xe[i] = strpos*sx0+2.*(sx0*strsig)*2.*(rand ()/(double)RAND_MAX - 0.5); else xe[i] = -strpos*sx0+2.*(sx0*strsig)*2.*(rand ()/(double)RAND_MAX - 0.5); ye[i] = rpipey*2.*(rand ()/(double)RAND_MAX - 0.5); } while (xe[i]*xe[i]/rpipex/rpipex + ye[i]*ye[i]/rpipey/rpipey >= 1); xpe[i] = 0.0; ype[i] = 0.0; } if (elflag==5) parabolic(); for (i=0; i= 1.0) ; xpr[i] = sx0 * u * sqrt(-2.0*log(s)/s); xpp[i] = sx0/betax0 * v * sqrt(-2.0*log(s)/s); do { u = 2.0 * rand () / (double)RAND_MAX - 1.0 ; v = 2.0 * rand () / (double)RAND_MAX - 1.0 ; s = u*u + v*v ; } while (s >= 1.0) ; ypr[i] = sy0 * u * sqrt(-2.0*log(s)/s); ypp[i] = sy0/betay0 * v * sqrt(-2.0*log(s)/s); do { u = 2.0 * rand () / (double)RAND_MAX - 1.0 ; v = 2.0 * rand () / (double)RAND_MAX - 1.0 ; s = u*u + v*v ; } while (s >= 1.0) ; sz[i] = z0_kick + sz0 * u * sqrt(-2.0*log(s)/s); dp[i] = dp0 * v * sqrt(-2.0*log(s)/s); if(isyn==3) sz[i] = sz0 * 4.0 * (rand () / (double)RAND_MAX - 0.5); xsum += xpr[i]; ysum += ypr[i]; xpsum += xpp[i]; ypsum += ypp[i]; } xave = xsum/(double)NPR; yave = ysum/(double)NPR; xpave = xpsum/(double)NPR; ypave = ypsum/(double)NPR; for (l=0; l= 1); xpe[i] = 0.0; ype[i] = 0.0; xeave += xe[i]; yeave += ye[i]; } if(elflag==3) for (i=0; i 0.5) xe[i] = strpos*sx0+2.*(sx0*strsig)*2.*(rand ()/(double)RAND_MAX - 0.5); else xe[i] = -strpos*sx0+2.*(sx0*strsig)*2.*(rand ()/(double)RAND_MAX - 0.5); ye[i] = rpipey*2.*(rand ()/(double)RAND_MAX - 0.5); } while (xe[i]*xe[i]/rpipex/rpipex + ye[i]*ye[i]/rpipey/rpipey >= 1); xpe[i] = 0.0; ype[i] = 0.0; xeave += xe[i]; yeave += ye[i]; } if(elflag==5) { parabolic(); for (i=0; i= 1); xpe[i+jmain] = 0.0; ype[i+jmain] = 0.0; xeave += xe[i+jmain]; yeave += ye[i+jmain]; } xeave /= (double)NEL; yeave /= (double)NEL; for (i=0; i= 0 && index < NBIN) { indexp[index][iind[index]] = i; xbin[index] += xpr[i]; xpbin[index] += xpp[i]; x2bin[index] += xpr[i]*xpr[i]; ybin[index] += ypr[i]; ypbin[index] += ypp[i]; y2bin[index] += ypr[i]*ypr[i]; ++iind[index]; ++ncount; } } for (i=0; iNPR1-1) { fprintf(stderr, "ERROR: too many particles in one slice !!"); exit(2); } if (iind[i] > 0) { xs[i] = xbin[i]/(double)iind[i]; ys[i] = ybin[i]/(double)iind[i]; xps[i] = xpbin[i]/(double)iind[i]; yps[i] = ypbin[i]/(double)iind[i]; npr[i] = (double)iind[i]/(double)NPR*npr0; if (iind[i] > 1) { sxs[i] = sqrt(x2bin[i]/(double)iind[i] - xs[i]*xs[i]); sys[i] = sqrt(y2bin[i]/(double)iind[i] - ys[i]*ys[i]); } else { sxs[i] = sx0_cd; sys[i] = sy0; } } else { xs[i] = 0.0; ys[i] = 0.0; xps[i] = 0.0; yps[i] = 0.0; npr[i] = 0.0; sxs[i] = sx0_cd; sys[i] = sy0; } } bin_frac = (double)ncount/(double)(NPR - npr_eff); } /****************************************************************************** *** Subroutine : wake_function *** *** Effect : Function that gives the wake function at a given *** *** location z < 0 (z is 'pos') *** *** Parameters : pos *** *** Gbl var used : *** *** Gbl var effect : *** *** Constants used : none *** *** Subrout. used : none *** ******************************************************************************/ double wake_func (double pos) { double field; field = (zetat*1e6*omegar*omegar/merit/omegabar)*exp(alpha*pos/C)*sin(omegabar*pos/C); return field; } /****************************************************************************** *** Subroutine : wake_function2 *** *** Effect : Function that gives the wake function at a given *** *** location z < 0 (z is 'pos') for a resistive wall *** *** Parameters : pos *** *** Gbl var used : *** *** Gbl var effect : *** *** Constants used : none *** *** Subrout. used : none *** ******************************************************************************/ double wake_reswall (double pos) { double sigma, field ; sigma = 5.4e13; //sigma = 5.4e17; /* copper conductivity in CGS [1/s] */ field = -1/PI/rpipex/rpipex/rpipex/2/PI/EPS0*sqrt(C/sigma)*1/sqrt(-pos)*circ; return field; } /****************************************************************************** *** Subroutine : longitudinal wake_function *** *** Effect : Function that gives the z-wake function at a given *** *** location z < 0 (z is 'pos') - single bunch in cavity *** *** Parameters : pos *** *** Gbl var used : *** *** Gbl var effect : *** *** Constants used : none *** *** Subrout. used : none *** ******************************************************************************/ double wake_funcz (double pos) { double field ; field = (rshz*1e6*2*alphaz)*exp(alphaz*pos/C)*(cos(omegabarz*pos/C) + alphaz/omegabarz*sin(omegabarz*pos/C)); return field; } /****************************************************************************** *** Subroutine : open_files *** *** Effect : Open files for storing calculation data. *** *** Parameters : none *** *** Gbl var used : *** *** Gbl var effect : *** *** Constants used : none *** *** Subrout. used : none *** ******************************************************************************/ open_files () { char bpr_filename[80] , /* bunch properties */ headtail_fname[80] , /* head tail distribution */ prphase_fname[80] , /* proton phase space */ pht_filename[80] , prot_filename[80] , samp_filename[80] ; strcpy (bpr_filename,filename) ; strcat (bpr_filename,"_prt.dat") ; bunch_pr = fopen(bpr_filename,"w"); strcpy (headtail_fname,filename) ; strcat (headtail_fname,"_hdtl.dat") ; headtail_pr = fopen(headtail_fname,"w"); strcpy (prphase_fname,filename) ; strcat (prphase_fname,"_prb.dat") ; prot_phase = fopen(prphase_fname,"w"); strcpy (pht_filename,filename) ; strcat (pht_filename,"_inph.dat") ; prot_inphase = fopen(pht_filename,"w"); strcpy (samp_filename,filename) ; strcat (samp_filename,"_sample.dat") ; samp_prot = fopen(samp_filename,"w"); strcpy (prot_filename,filename) ; strcat (prot_filename,"_bunchds.dat") ; prot_distr = fopen(prot_filename,"w"); } /****************************************************************************** **** PIC module **** ******************************************************************************/ void* xmalloc (size_t size) { void* tmp; tmp=malloc(size); if (tmp!=NULL) { return tmp; } fprintf(stderr,"ERROR: not enough memory\n"); exit(1); } /* Distributes the beam particles for field calculation */ void particles_distribute(GRID *grid,double x[],double y[],int n, int i_which, double n_macro) { int i,j,i1,i2,n_h,j1,j2,j3,j4; double delta_x_i,delta_y_i,offset_x,offset_y; double h_x,h_y; double *rho,tmp; int n_x,n_y; int distribute_in=0,distribute_out=0; n_x=grid->n_cell_x; n_y=grid->n_cell_y; if (i_which==1) { rho=grid->rho1; } else { rho=grid->rho2; } delta_x_i=1.0/grid->delta_x; delta_y_i=1.0/grid->delta_y; offset_x=grid->offset_x; offset_y=grid->offset_y; for (i=0;idistribution) { case 1: for (i=0;i=0)&&(i1=0)&&(i2=0)&&(i1=0)&&(i2= i__2 : i2 <= i__2; i2 += i__3) { if (i2 < i2rev) { i__4 = i2 + ip1 - 2; for (i1 = i2; i1 <= i__4; i1 += 2) { i__5 = ip3; i__6 = ip2; for (i3 = i1; i__6 < 0 ? i3 >= i__5 : i3 <= i__5; i3 += i__6) { i3rev = i2rev + i3 - i2; tempr = data[i3]; tempi = data[i3 + 1]; data[i3] = data[i3rev]; data[i3 + 1] = data[i3rev + 1]; data[i3rev] = tempr; data[i3rev + 1] = tempi; } } } ibit = ip2 / 2; L1: if (ibit >= ip1 && i2rev > ibit) { i2rev -= ibit; ibit /= 2; goto L1; } i2rev += ibit; } ifp1 = ip1; L2: if (ifp1 < ip2) { ifp2 = ifp1 << 1; theta = isign * 6.28318530717959 / (ifp2 / ip1); /* Computing 2nd power */ d__1 = sin(theta * .5); wpr = d__1 * d__1 * -2.; wpi = sin(theta); wr = 1.; wi = 0.; i__3 = ifp1; i__2 = ip1; for (i3 = 1; i__2 < 0 ? i3 >= i__3 : i3 <= i__3; i3 += i__2) { i__4 = i3 + ip1 - 2; for (i1 = i3; i1 <= i__4; i1 += 2) { i__6 = ip3; i__5 = ifp2; for (i2 = i1; i__5 < 0 ? i2 >= i__6 : i2 <= i__6; i2 += i__5) { k1 = i2; k2 = k1 + ifp1; tempr = wr * data[k2] - wi * data[k2 + 1]; tempi = wr * data[k2 + 1] + wi * data[k2]; data[k2] = data[k1] - tempr; data[k2 + 1] = data[k1 + 1] - tempi; data[k1] += tempr; data[k1 + 1] += tempi; } } wtemp = wr; wr = wr * wpr - wi * wpi + wr; wi = wi * wpr + wtemp * wpi + wi; } ifp1 = ifp2; goto L2; } nprev = n * nprev; } return; } /* This routine calculates the potentials with the help of the fast fourier transform */ void fold_fft(double *rho1,double *rho2,double *dist,double *phi1, double *phi2,double temp[],int n_x,int n_y) { FILE *datei; double help; int i1,i2,i,j,j0; int nn[2]; double eps=1e-5; /* return no field */ /* fill array with charge */ for (i=0;i0)&&(anorm0)&&(anormn_cell_x; n_y=grid->n_cell_y; grid->min_x=-((float)n_x-2)/((float)n_x)*cut_x; grid->max_x=((float)n_x-2)/((float)n_x)*cut_x; grid->min_y=-((float)n_y-2)/((float)n_y)*cut_y; grid->max_y=((float)n_y-2)/((float)n_y)*cut_y; grid->offset_x=0.5*(float)n_x; grid->offset_y=0.5*(float)n_y; grid->delta_x=2.0*cut_x/((double)n_x); grid->delta_y=2.0*cut_y/((double)n_y); factor=1.0/(grid->delta_x*grid->delta_y); switch (grid->integration_method) { case 1: case 3: phi0=factor* f_potential_2(0.0,0.0,0.5*grid->delta_x,0.5*grid->delta_y); for (i1=0;i1n_cell_x;i1++) { for (i2=0;i2n_cell_y;i2++) { j=i1*grid->n_cell_y+i2; x0=(double)i1; y0=(double)i2; grid->dist[j]=factor* f_potential_2(x0*grid->delta_x,y0*grid->delta_y, 0.5*grid->delta_x,0.5*grid->delta_y) -phi0; } } break; case 2: phi0=factor* f_potential_2(0.0,0.0,0.5*grid->delta_x,0.5*grid->delta_y); for (i1=0;i1dist[i1]=0.0; } for (i1=0;i1dist[j0]=factor* f_potential_2(x0*grid->delta_x,y0*grid->delta_y, 0.5*grid->delta_x,0.5*grid->delta_y) -phi0; if (i2>=1) { grid->dist[2*((i1+1)*2*n_y-i2)]=grid->dist[j0]; if (i1>=1) { grid->dist[2*(2*n_y*(2*n_x-i1)+i2)]=grid->dist[j0]; grid->dist[2*(2*n_y*(2*n_x-i1)+2*n_y-i2)]=grid->dist[j0]; } } else { if (i1>=1) { grid->dist[2*(2*n_y*(2*n_x-i1)+i2)]=grid->dist[j0]; } } } } nn[0]=2*n_y; nn[1]=2*n_x; i1=2; i2=1; //fourtrans(grid->dist,nn,i2); fourtrans2(grid->dist,nn,i2); break; } } GRID* grid_init_comp (int n_x,int n_y,int integration_method) { GRID *grid; int i1,i2,j; grid=(GRID*)xmalloc(sizeof(GRID)); grid->n_cell_x=n_x; grid->n_cell_y=n_y; grid->integration_method=integration_method; grid->distribution=2; grid->interpolation=2; grid->rho1=(double*)xmalloc(sizeof(double)*n_x*n_y); grid->rho2=(double*)xmalloc(sizeof(double)*n_x*n_y); grid->phi1=(double*)xmalloc(sizeof(double)*n_x*n_y); grid->phi2=(double*)xmalloc(sizeof(double)*n_x*n_y); if (integration_method==2) { grid->dist=(double*)xmalloc(sizeof(double)*n_x*n_y*8); grid->temp=(double*)xmalloc(sizeof(double)*n_x*n_y*8); } else { grid->dist=(double*)xmalloc(sizeof(double)*n_x*n_y); } return grid; } /* This routine folds the charge distribution and the greensfunction of a grid to get the potentials. */ void foldfields (double *rho,double *dist,double *phi,int n_x,int n_y) { double s; int i1,i2,i3,i4,j; int field_typ=1; float eps=1e-5; for (i1=0;i1integration_method) { case 1: foldfields(grid->rho1,grid->dist,grid->phi1,grid->n_cell_x,grid->n_cell_y); foldfields(grid->rho2,grid->dist,grid->phi2,grid->n_cell_x,grid->n_cell_y); break; case 2: fold_fft(grid->rho1,grid->rho2,grid->dist,grid->phi1,grid->phi2, grid->temp,grid->n_cell_x,grid->n_cell_y); break; } } void particles_move(GRID *grid,double energy,double ics0[], double yps0[],double icsp[],double ypsp[], int nn,int i_which,double step,double scale) { static float esum=0.0; static int sum_photon=0; static long int n_low=0; int i,i1,i2,l,k; double ax,ay; double phi1_x,phi2_x,phi3_x,phi1_y,phi2_y,phi3_y,h_x,h_y; int n_x,n_y; register int j; register double h,h_p; double *phi; double delta_x_i,delta_y_i,offset_x,offset_y,delta_x,delta_y; double xx,yy,xxp,yyp; double min_x,max_x,min_y,max_y,e_tmp; double k11,k12,k13,k14,k21,k22,k23,k24,k31,k32,k33,k34,k41,k42,k43,k44; double velprx, velpry, bux, buy, bmod; double p; min_x=grid->min_x; max_x=grid->max_x; min_y=grid->min_y; max_y=grid->max_y; delta_x=grid->delta_x; delta_y=grid->delta_y; delta_x_i=1.0/grid->delta_x; delta_y_i=1.0/grid->delta_y; offset_x=grid->offset_x; offset_y=grid->offset_y; n_x=grid->n_cell_x; n_y=grid->n_cell_y; if (i_which==1){ phi=grid->phi2; } else{ phi=grid->phi1; } switch (grid->interpolation){ case 1: for (i=0;i0)&&(i20)){ j=i1*n_y+i2; ax = (phi[j+n_y] -phi[j-n_y])*step /(2.0*delta_x*energy); ay = (phi[j+1] -phi[j-1])*step /(2.0*delta_y*energy); icsp[i] -= ax*scale; ypsp[i] -= ay*scale; } else { fprintf(stderr,"WARNING: particle is not on the grid\n"); } ics0[i] += icsp[i]*step; yps0[i] += ypsp[i]*step; } break; case 2: for (i=0;i=rpipex || fabs(yps0[i])>=rpipey) { do { if(i_dip==1) icsp[i] = 0.0; else icsp[i] = 1.3e6*rand () / (double)RAND_MAX ; ypsp[i] = 1.3e6*rand () / (double)RAND_MAX ; } while (sqrt(icsp[i]*icsp[i]+ypsp[i]*ypsp[i]) >= 1.3e6) ; if(fabs(ics0[i])>=rpipex) ics0[i] = fabs(ics0[i])/ics0[i]*rpipex; if(fabs(yps0[i])>=rpipey) yps0[i] = fabs(yps0[i])/yps0[i]*rpipey; icsp[i] *= -fabs(ics0[i])/ics0[i]; ypsp[i] *= -fabs(yps0[i])/yps0[i]; } } if(elflag==2) { if (ics0[i]*ics0[i]/rpipex/rpipex+yps0[i]*yps0[i]/rpipey/rpipey>=1.) { do { if(i_dip==1) icsp[i] = 0.0; else icsp[i] = 1.3e6*rand () / (double)RAND_MAX ; ypsp[i] = 1.3e6*rand () / (double)RAND_MAX ; } while (sqrt(icsp[i]*icsp[i]+ypsp[i]*ypsp[i]) >= 1.3e6) ; ics0[i] = rpipex*cos(atan2(yps0[i],ics0[i])); yps0[i] = rpipey*sin(atan2(yps0[i],ics0[i])); icsp[i] *= -fabs(ics0[i])/ics0[i]; ypsp[i] *= -fabs(yps0[i])/yps0[i]; } } } xx=ics0[i]; yy=yps0[i]; xxp=icsp[i]; yyp=ypsp[i]; if ((xx>min_x)&&(xxmin_y)&&(yy Reading data from configuration file.\n") ; read_data () ; printf ("--> Init values.\n") ; init_values () ; printf ("\t\t Complete\n") ; open_files (); printf ("--> Generation of e- cloud and bunch distributions.\n"); initialization (); printf ("\t\t Complete\n") ; if(i_method==2) { printf ("--> Generation of the grid for the particle-in-cell.\n"); printf (" calculation of the fields \n"); grid=grid_init_comp(128,128,2); if(i_dip==2){ gridextx = 3.*rpipex; gridexty = 3.*rpipey; } else{ gridextx = 1.1*rpipex; gridexty = 1.1*rpipey; } grid_init_phys(grid,gridextx,gridexty); printf ("\t\t Complete\n") ; } fprintf (prot_inphase, "electrons = %d protons = %d bunch slices = %d\n",NEL,NPR,NBIN); fprintf (prot_inphase, "electron sizes -> sx = %g sy = %g\n",sx2,sy2); fprintf (prot_inphase, "beam sizes -> sx = %g sy = %g\n",sx0,sy0); fprintf (prot_inphase, "e- cloud sizes -> rpipex = %g rpipey = %g\n",rpipex,rpipey); fprintf (prot_inphase, "maximum tune shifts -> dt_x = %g dt_y = %g\n",dtunxmax,dtunymax); ili = 0; /* here the loop over the total number of steps (kicks * turns) starts... */ for(it=1; it<=nstep; it++) { if(npr_eff==NPR) exit(2); xave = 0.0; yave = 0.0; zave = 0.0; xpave = 0.0; ypave = 0.0; dpmean = 0.0; sxsize = 0.0; sysize = 0.0; szsize = 0.0; xemit = 0.0; yemit = 0.0; corryz = 0.0; invary = 0.0; invarx = 0.0; for(imain=npr_eff; imain0.1) for(jmain=npr_eff; jmain=0; jmain--) { sx = sxs[jmain]; sy = sys[jmain]; xoff = xs[jmain]; yoff = ys[jmain]; npr1 = npr[jmain]; /* offsets and rms sizes of the slices are printed to file */ if(it % n_diag == 0 || it == 1) fprintf(headtail_pr, "%13.8e %13.8e %13.8e %13.8e %13.8e\n",zs[jmain],npr1*xoff,npr1*yoff,npr1*sx,npr1*sy); /* here the loop over the beam-kicked electrons starts... if the soft-Gaussian model is applied */ if(i_dip==2) new_fraction (); switch (i_method) { case 1: for(imain=0; imain0) { for(kmain=0; kmain0 && iwake_flag == 1) { for (kmain=0; kmainjmain; mmain--) { zi = (jmain-mmain)*zstep; switch(i_pipe) { case 0: kick_x += wakefac*npr[mmain]*xs[mmain]*wake_func(zi); kick_y += wakefac*npr[mmain]*ys[mmain]*wake_func(zi); kick_z += wakefac*npr[mmain]*wake_funcz(zi); break; case 1: kick_x += wakefac*npr[mmain]*0.41/1.23*(xs[mmain]-xpr[ind])*wake_func(zi); kick_y += wakefac*npr[mmain]*0.82/1.23*(ys[mmain]+0.5*ypr[ind])*wake_func(zi); kick_z += wakefac*npr[mmain]*wake_funcz(zi); case 2: kick_x += 0.; kick_y += wakefac*npr[mmain]*ys[mmain]*wake_func(zi); kick_z += wakefac*npr[mmain]*wake_funcz(zi); case 3: kick_x += wakefac*npr[mmain]*xs[mmain]*wake_reswall(zi); kick_y += wakefac*npr[mmain]*ys[mmain]*wake_reswall(zi); kick_z += wakefac*npr[mmain]*wake_funcz(zi); break; } //if(it==1 && kmain==0 && mmain==NBIN-1) //fprintf(prot_distr,"%13.8e %13.8e\n",zi,wake_func(merit,zetat,zi)); } } else kick_z = -9.57e-10 + 8.e-11*2.*(rand()/(double)RAND_MAX - 0.5); xpp[ind] += kick_x; ypp[ind] += kick_y; dp[ind] += kick_z; } } if (i_space == 1) /* { for(kmain=0; kmain0.1) for(jmain=npr_eff; jmain= 0 && sl_ind < NBIN) { cent_x = xs[sl_ind]; cent_y = ys[sl_ind]; cent_xp = xps[sl_ind]; cent_yp = yps[sl_ind]; dtunx = dtunxmax*sx0_cd*sx0_cd/sxs[sl_ind]/(sxs[sl_ind]+sys[sl_ind])*npr[sl_ind]/npr0*sqrt(2*PI)*sz0/zstep; dtuny = dtunymax*sy0*sy0/sys[sl_ind]/(sxs[sl_ind]+sys[sl_ind])*npr[sl_ind]/npr0*sqrt(2*PI)*sz0/zstep; //printf("dtunx = %lf", dtunx); //printf("dtuny = %lf", dtuny); //dtunx = dtunxmax*sx0_cd*sx0_cd/sxs[sl_ind]/(sxs[sl_ind]+sys[sl_ind])*lambda; //dtuny = dtunymax*sy0*sy0/sys[sl_ind]/(sxs[sl_ind]+sys[sl_ind])*lambda; } else { if (sl_ind < 0) { cent_x = xs[0]; cent_y = ys[0]; cent_xp = xps[0]; cent_yp = yps[0]; dtunx = dtunxmax*sx0_cd*sx0_cd/sxs[0]/(sxs[0]+sys[0])*npr[0]/npr0*sqrt(2*PI)*sz0/zstep; dtuny = dtunymax*sy0*sy0/sys[0]/(sxs[0]+sys[0])*npr[0]/npr0*sqrt(2*PI)*sz0/zstep; //dtunx = dtunxmax*sx0_cd*sx0_cd/sxs[0]/(sxs[0]+sys[0])*lambda; //dtuny = dtunymax*sy0*sy0/sys[0]/(sxs[0]+sys[0])*lambda; } else { cent_x = xs[NBIN-1]; cent_y = ys[NBIN-1]; cent_xp = xps[NBIN-1]; cent_yp = yps[NBIN-1]; dtunx = dtunxmax*sx0_cd*sx0_cd/sxs[NBIN-1]/(sxs[NBIN-1]+sys[NBIN-1])*npr[NBIN-1]/npr0*sqrt(2*PI)*sz0/zstep; dtuny = dtunymax*sy0*sy0/sys[NBIN-1]/(sxs[NBIN-1]+sys[NBIN-1])*npr[NBIN-1]/npr0*sqrt(2*PI)*sz0/zstep; //dtunx = dtunxmax*sx0_cd*sx0_cd/sxs[NBIN-1]/(sxs[NBIN-1]+sys[NBIN-1])*lambda; //dtuny = dtunymax*sy0*sy0/sys[NBIN-1]/(sxs[NBIN-1]+sys[NBIN-1])*lambda; } } xx1 = cos(2*PI*dtunx/(double)nkick)*(x0 - cent_x) - betax0*sin(2*PI*dtunx/(double)nkick)*(xp0 - cent_xp) + cent_x; xxp1 = 1/betax0*sin(2*PI*dtunx/(double)nkick)*(x0 - cent_x) + cos(2*PI*dtunx/(double)nkick)*(xp0 - cent_xp) + cent_xp; yy1 = cos(2*PI*dtuny/(double)nkick)*(ips0 - cent_y) - betay0*sin(2*PI*dtuny/(double)nkick)*(yp0 - cent_yp) + cent_y; yyp1 = 1/betay0*sin(2*PI*dtuny/(double)nkick)*(ips0 - cent_y) + cos(2*PI*dtuny/(double)nkick)*(yp0 - cent_yp) + cent_yp; } if(i_rot==1) { xx1 = cos(2*PI*dtunx/(double)nkick)*(x0 - cent_x) - betax0*sin(2*PI*dtunx/(double)nkick)*(xp0 - cent_xp) + cent_x; xxp1 = 1/betax0*sin(2*PI*dtunx/(double)nkick)*(x0 - cent_x) + cos(2*PI*dtunx/(double)nkick)*(xp0 - cent_xp) + cent_xp; yy1 = cos(2*PI*dtuny/(double)nkick)*(ips0 - cent_y) - betay0*sin(2*PI*dtuny/(double)nkick)*(yp0 - cent_yp) + cent_y; yyp1 = 1/betay0*sin(2*PI*dtuny/(double)nkick)*(ips0 - cent_y) + cos(2*PI*dtuny/(double)nkick)*(yp0 - cent_yp) + cent_yp; } break; case 2: lambda = 2*(rand ()/(double)RAND_MAX - 0.5); xx1 = cos(2*PI*dtunxmax*lambda/(double)nkick)*x0 - betax0*sin(2*PI*dtunxmax*lambda/(double)nkick)*xp0; xxp1 = 1/betax0*sin(2*PI*dtunxmax*lambda/(double)nkick)*x0 + cos(2*PI*dtunxmax*lambda/(double)nkick)*xp0; yy1 = cos(2*PI*dtunymax*lambda/(double)nkick)*ips0 - betay0*sin(2*PI*dtunymax*lambda/(double)nkick)*yp0; yyp1 = 1/betay0*sin(2*PI*dtunymax*lambda/(double)nkick)*ips0 + cos(2*PI*dtunymax*lambda/(double)nkick)*yp0; break; } if(i_oct==1) { epsx = (x0*x0 + xp0*xp0*betax0*betax0)/betax0; epsy = (ips0*ips0 + yp0*yp0*betay0*betay0)/betay0; /* Remember that epsx = 2*actx; epsy = 2*acty */ cxa = cos(2*PI*(htune + qpxx*epsx + qpxy*epsy)/(double)nkick); sxa = sin(2*PI*(htune + qpxx*epsx + qpxy*epsy)/(double)nkick); cya = cos(2*PI*(vtune + qpxy*epsx + qpyy*epsy)/(double)nkick); sya = sin(2*PI*(vtune + qpxy*epsx + qpyy*epsy)/(double)nkick); } r11 = cxa*cos(2*PI*csix*dpp0/(double)nkick) - sxa * sin(2*PI*csix*dpp0/(double)nkick); r12 = betax0*(sxa*cos(2*PI*csix*dpp0/(double)nkick) + cxa * sin(2*PI*csix*dpp0/(double)nkick)); r21 = - r12/betax0/betax0; r22 = r11; /* horizontal coordinates transformation - linear and with chromatic rotation */ if(i_incoh==0) { if(ishift!=1) { sl_ind = (long)floor((szz0 + 2.*szsize)/(4.*szsize)*NBIN); if (sl_ind >= 0 && sl_ind < NBIN) { cent_x = xs[sl_ind]; cent_y = ys[sl_ind]; cent_xp = xps[sl_ind]; cent_yp = yps[sl_ind]; } else { if (sl_ind < 0) { cent_x = xs[0]; cent_y = ys[0]; cent_xp = xps[0]; cent_yp = yps[0]; } else { cent_x = xs[NBIN-1]; cent_y = ys[NBIN-1]; cent_xp = xps[NBIN-1]; cent_yp = yps[NBIN-1]; } } } xx1 += -cent_x; xxp1 += -cent_xp; yy1 += -cent_y; yyp1 += -cent_yp; } /* next line introduces the possibility to have linear coupling in the machine */ if(i_coupl==1) { xxp1 -= key_xy*yy1; yyp1 -= key_xy*xx1; } xpr[kmain] = xxp1 * r12 + xx1 * r11; xpp[kmain] = xxp1 * r22 + xx1 * r21; r33 = cya*cos(2*PI*csiy*dpp0/(double)nkick) - sya * sin(2*PI*csiy*dpp0/(double)nkick); r34 = betay0*(sya*cos(2*PI*csiy*dpp0/(double)nkick) + cya * sin(2*PI*csiy*dpp0/(double)nkick)); r43 = - r34/betay0/betay0; r44 = r33; /* vertical coordinates transformation - linear and with chromatic rotation */ ypr[kmain] = yyp1 * r34 + yy1 * r33; ypp[kmain] = yyp1 * r44 + yy1 * r43; /* if one chooses to simulate synchrotron motion too, the longitudinal coordinates are transformed according to a linear focusing law */ switch(isyn){ case 0: break; case 1: sz[kmain] = csa * szz0 - eta * C/omegas * ssa * dpp0; dp[kmain] = csa * dpp0 + omegas/eta/C * ssa * szz0; break; case 2: sz[kmain] = szz0 - eta*circ*dpp0; break; case 3: sz[kmain] = szz0 - eta*circ*dpp0; if (fabs(sz[kmain]) >= 2*sz0) { sz[kmain] = 4*sz0*fabs(sz[kmain])/sz[kmain] - sz[kmain]; dp[kmain] = -dp[kmain]; } break; case 4: kp1 = coeff3*sin(coeff4*szz0)*tstepp; k1 = coeff5*dpp0*tstepp; kp2 = coeff3*sin(coeff4*(szz0 + k1/2))*tstepp; k2 = coeff5*(dpp0 + kp1/2)*tstepp; kp3 = coeff3*sin(coeff4*(szz0 + k2/2))*tstepp; k3 = coeff5*(dpp0 + kp2/2)*tstepp; kp4 = coeff3*sin(coeff4*(szz0 + k3))*tstepp; k4 = coeff5*(dpp0 + kp3)*tstepp; sz[kmain] = szz0 + k1/6 + k2/3 + k3/3 + k4/6; dp[kmain] = dpp0 + kp1/6 + kp2/3 + kp3/3 + kp4/6; break; case 5: kp1 = (coeff3*sin(coeff4*szz0) + coeffspch*(szz0 - zave)/(szsize*szsize*szsize)*exp(-(szz0 - zave)*(szz0 - zave)/(2*szsize*szsize)))*tstepp; k1 = coeff5*dpp0*tstepp; kp2 = (coeff3*sin(coeff4*(szz0 + k1/2)) + coeffspch*(szz0 - zave)/(szsize*szsize*szsize)*exp(-(szz0 - zave)*(szz0 - zave)/(2*szsize*szsize)))*tstepp; k2 = coeff5*(dpp0 + kp1/2)*tstepp; kp3 = (coeff3*sin(coeff4*(szz0 + k2/2)) + coeffspch*(szz0 - zave)/(szsize*szsize*szsize)*exp(-(szz0 - zave)*(szz0 - zave)/(2*szsize*szsize)))*tstepp; k3 = coeff5*(dpp0 + kp2/2)*tstepp; kp4 = (coeff3*sin(coeff4*(szz0 + k3)) + coeffspch*(szz0 - zave)/(szsize*szsize*szsize)*exp(-(szz0 - zave)*(szz0 - zave)/(2*szsize*szsize)))*tstepp; k4 = coeff5*(dpp0 + kp3)*tstepp; sz[kmain] = szz0 + k1/6 + k2/3 + k3/3 + k4/6; dp[kmain] = dpp0 + kp1/6 + kp2/3 + kp3/3 + kp4/6; break; } if (kmain % 300==0 && it % n_diag==0) fprintf(prot_phase,"%d \t%13.8e\t%13.8e\t%13.8e\t%13.8e\t%13.8e\t%13.8e\n",kmain,xpr[kmain],xpp[kmain],ypr[kmain],ypp[kmain],sz[kmain],dp[kmain]); else if (kmain==NPR-1) fprintf(samp_prot,"%d \t%13.8e\t%13.8e\t%13.8e\t%13.8e\t%13.8e\t%13.8e\n",it,xpr[kmain],xpp[kmain],ypr[kmain],ypp[kmain],sz[kmain],dp[kmain]); if(fabs(xpr[kmain])>=rpipex || fabs(ypr[kmain])>=rpipey) { xpr[kmain] = xpr[npr_eff]; ypr[kmain] = ypr[npr_eff]; sz[kmain] = sz[npr_eff]; xpp[kmain] = xpp[npr_eff]; ypp[kmain] = ypp[npr_eff]; dp[kmain] = dp[npr_eff]; ++npr_eff; } } if (it % n_diag == 0) { if(isyn!=3) verteil (sz,NPR-npr_eff,hdp,NDC,indxp); else verteil2 (sz,NPR-npr_eff,hdp,NDC,indxp); for (lmain=0; lmain %g \n",(double)(npr_eff)/(double)(NPR)*100); fclose(bunch_pr); fclose(headtail_pr); fclose(prot_phase); fclose(prot_inphase); fclose(prot_distr); fclose(samp_prot); }