1 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
---|
2 | c WARNING c |
---|
3 | c Warning for change of System and Operating system we have to do c |
---|
4 | c one change of the variable kvers in main program ligne 282 c |
---|
5 | c utilised in ddate c |
---|
6 | c----------------------------------------------------------------------c |
---|
7 | c c |
---|
8 | c VERSION V5.03 c |
---|
9 | c c |
---|
10 | c October 90 Orsay V4.00 c |
---|
11 | c Release Mai 92 V4.10 c |
---|
12 | c Release June 92 V4.20 c |
---|
13 | c Release October 92 V4.21 c |
---|
14 | c Release December 92 V4.22 c |
---|
15 | c Release July 93 V4.30 c |
---|
16 | c Release March 94 V4.31 c |
---|
17 | c Release April 94 V4.32 c |
---|
18 | c Release June 95 V5.00 c |
---|
19 | c Release June 96 V5.01 c |
---|
20 | c Release July 96 V5.02 VMS / UNIX c |
---|
21 | c Release January 97 V5.03 VMS / UNIX c |
---|
22 | c Release September 2000 V5.03 Multi plateforme sans cernlib c |
---|
23 | c Release April 2010 V5.03 devient version officielle LAL c |
---|
24 | c c |
---|
25 | c c |
---|
26 | c B. Mouton LAL Bat 200 91405 Orsay France c |
---|
27 | c c |
---|
28 | c MOUTON@LAL.IN2P3.FR c |
---|
29 | c c |
---|
30 | c Version ALPHA c |
---|
31 | c This version of PARMELA is a pachwork of 3 versions c |
---|
32 | c one of from Lloyd Young (version of July 88) of Los Alamos, c |
---|
33 | c another of KT Mc Donald of Brookhaven and one of SLAC with c |
---|
34 | c some local modifications, with modifications and news c |
---|
35 | c routines from H Liu of IHEP Beijing,with news routines from c |
---|
36 | c CEBAF. c |
---|
37 | c c |
---|
38 | c---PhAse and Radial Motion in Electron Linear Accelerators c |
---|
39 | c---(and associated transport systems) c |
---|
40 | c c |
---|
41 | c on VAX VMS you have to increase your Paging file quota to 20000 c |
---|
42 | c c |
---|
43 | c $for/noop parmelav503 c |
---|
44 | c $link/exe=parmelav5 parmvelav503,'lib$' c |
---|
45 | c c |
---|
46 | c # Makefile for PARMELA cas un seul fichier c |
---|
47 | c # on Alpha OSF1 c |
---|
48 | c # c |
---|
49 | c parmela : parmelav5.f c |
---|
50 | c f77 -g -O0 parmelav5.f -o parmelav5 ${CERNLIBS} c |
---|
51 | c c |
---|
52 | c # on HP c |
---|
53 | c parmela : parmelav5.f c |
---|
54 | c f77 +ppu -K -g -w +E1 -O parmelav503.f -o parmelav5hp ${CERNLIBS}c |
---|
55 | c c |
---|
56 | c----------------------------------------------------------------------c |
---|
57 | c c |
---|
58 | c Pour le cas de la version eclatee 07/2000 c |
---|
59 | c #!//bin/ksh c |
---|
60 | c # gmake sans cernlib c |
---|
61 | c programe="parmela" c |
---|
62 | c c |
---|
63 | c for files in *.f c |
---|
64 | c do c |
---|
65 | c objname="${files:%.f}" c |
---|
66 | c echo $objname c |
---|
67 | c objname=${objname}.o c |
---|
68 | c liste=${liste}" ${files}" c |
---|
69 | c #echo $liste c |
---|
70 | c echo ${objname}:${files} >>temp1 c |
---|
71 | c echo "\t \$(FF) -c \$(FFLAGS) "${files} >>temp1 c |
---|
72 | c echo " ">>temp1 c |
---|
73 | c done c |
---|
74 | c c |
---|
75 | c echo "${programe}: \$(OBJ)" >>temp1 c |
---|
76 | c #echo "\t \$(FF) \$(FFLAGS) \$(OBJ) -o "${programe}"\$(VER)" c |
---|
77 | c \${CERNLIBS} >>temp1 c |
---|
78 | c echo "\t \$(FF) \$(FFLAGS) \$(OBJ) -o "${programe}"\$(VER)" c |
---|
79 | c >>temp1 c |
---|
80 | c echo " " >>temp1 c |
---|
81 | c echo "clean:" >>temp1 c |
---|
82 | c echo "\t rm *.o">>temp1 c |
---|
83 | c c |
---|
84 | c # ecriture dans le fichier de makefile definitif c |
---|
85 | c c |
---|
86 | c ## Pour HP c |
---|
87 | c echo FF=f77 >>makefile c |
---|
88 | c #echo FFLAGS=+ppu +U77 -K -g -w +E1 -O >>makefile c |
---|
89 | c echo FFLAGS= -K -w +E1 +Oall >>makefile c |
---|
90 | c echo " " >>makefile c |
---|
91 | c echo "VER=alonehpopt" >>makefile c |
---|
92 | c c |
---|
93 | c ## Pour OSF c |
---|
94 | c #echo FF=f77 >>makefile c |
---|
95 | c #echo FFLAGS= -g -O2 >>makefile c |
---|
96 | c c |
---|
97 | c echo " " >>makefile c |
---|
98 | c echo SRC=${liste} >>makefile c |
---|
99 | c echo " " >>makefile c |
---|
100 | c echo "OBJ=\$(SRC:.f=.o)" >>makefile c |
---|
101 | c echo " " >>makefile c |
---|
102 | c cat temp1 >> makefile c |
---|
103 | c rm -f temp1 c |
---|
104 | c #c'est fini c |
---|
105 | c----------------------------------------------------------------------c |
---|
106 | c c |
---|
107 | c modification des call scheff pour input0 02/92 c |
---|
108 | c modification dans emit90 pour kk=0 05/92 --> version 4.10 c |
---|
109 | c modification dans elips pour negative square root 05/92 c |
---|
110 | c outlal, emit90 and elips are in double precision 05/92 c |
---|
111 | c out6 is in double precision 06/92 c |
---|
112 | c version 4.20 02/06/92 c |
---|
113 | c modification dans out6dp boucle 9996 remplace cor par cord c |
---|
114 | c version 4.21 c |
---|
115 | c version 4.22 ip.eq.999 for imp.dat cell and trwave c |
---|
116 | c modification dans cellfeld rcir en rcur 04/93 c |
---|
117 | c version 4.30 suppression of calls hbook c |
---|
118 | c version 4.31 24/03/94 c |
---|
119 | c pvol is in double precision pvoldp 03/94 c |
---|
120 | c pvoldp supressed icall 03/94 c |
---|
121 | c intro de cord(6,1)=cor(6,1) in main after 125 continue 03/94c |
---|
122 | c scheff2 loop 20 i=2,nr ----> i=2,nr1 03/94 c |
---|
123 | c in out6dp change in format 220 zmax = f8.4 in c |
---|
124 | c zmax = f8.2 03/94 c |
---|
125 | c version 4.32 18/04/94 c |
---|
126 | c introduction scaling factor on foclal card 04/94 c |
---|
127 | c modif /pcord/ intro nupar 05/94 c |
---|
128 | c nspl=10000 ----> 10001 11/94 c |
---|
129 | c nsplc=1497 ----> 30005=10001*3+2 11/94 c |
---|
130 | c in zoutlal replaced in data bcd : design by trwave 11/94 c |
---|
131 | c renumber subroutine input 27/03/95 c |
---|
132 | c input 11 27/03/95 c |
---|
133 | c histogram1 and histogram2 28/03/95 c |
---|
134 | c renumber main program 29/03/95 c |
---|
135 | c subroutine poisson 30/03/95 c |
---|
136 | c input 12 31/03/95 c |
---|
137 | c subroutine sextupole 11/04/95 c |
---|
138 | c modification of space charge card 11/04/95 c |
---|
139 | c new parameter iprog ; c |
---|
140 | c scheff beami iprog point and subroutine dependent c |
---|
141 | c modification extension infld in sffld 01/06/95 c |
---|
142 | c zout introduction de angle pour signe arctangente 04/09/95 c |
---|
143 | c modification input11 write ngood et tmax=2.5 sigma 12/09/95 c |
---|
144 | c modif dans input11 du calcul de degmax c |
---|
145 | c remplace tmax par tfac 28/09/95 c |
---|
146 | c intro dans parmela routine gpindp plus dans cernlib 28/09/95c |
---|
147 | c deplacement de l'open de parmtrch cas ip=3 19/02/96 c |
---|
148 | c modification input 4 26/06/96 c |
---|
149 | c - 8 parameters one first card input 4 c |
---|
150 | c modification input11 05/07/96 two news parameters angle of c |
---|
151 | c injection of laser light and distane focal point cathode c |
---|
152 | c modification input10 23/07/96 two news parameters angle of c |
---|
153 | c injection of laser light and distane focal point cathode c |
---|
154 | c 13/12/96 modif zout pour trwave (phase) c |
---|
155 | c 13/12/96 modif fieldlal pour nchp.lt.nspl c |
---|
156 | c 13/01/97 modif fieldlal pour free format c |
---|
157 | c 13/01/97 modif common/image/ c |
---|
158 | c 25/02/97 intro dans subroutine input pour input15 c |
---|
159 | c dppi,dpepsout double precision et modif de l'appel gpindp c |
---|
160 | c 03/09/97 mise sous Unix des statistiques elapsed cpu time c |
---|
161 | c 04/09/97 version pour HP (c_hp) c |
---|
162 | c 14/11/97 nettoyage des variables non utilisees c |
---|
163 | c intro du call date_and_time sous VMS pour an 2000 c |
---|
164 | c tous les open avec status='scratch' sans filename c |
---|
165 | c 03/02/98 modifation de zout et zoutlal pour cas solenoid c |
---|
166 | c 18/09/98 introduction de iadp a la place ngood dans out6dp c |
---|
167 | c write(ndes2), write(ndes1) c |
---|
168 | c 29/03/00 reprise version eclate a partir du source non c |
---|
169 | c non eclate du 31/05/99 pour correction d'un bug c |
---|
170 | c dans une version ephemere dite v503b c |
---|
171 | c 18/07/00 abandon de la cernlib c |
---|
172 | c suppression de call histo dans main et input c |
---|
173 | c 20/07/00 Modification de ddate.f pour essayer une methode c |
---|
174 | c universelle de date, modification du main c |
---|
175 | c introduction de kvers, introduction de la c |
---|
176 | c subroutine datime.f machine dependante c |
---|
177 | c suppression de l'utilisation de la cernlib c |
---|
178 | c 13/09/00 Introduction des include 'filename.h' c |
---|
179 | c Les fichiers .h contiennent les parameters et les c |
---|
180 | c common. c |
---|
181 | c Modification des noms de variable dans scwall c |
---|
182 | c Modification npoinz en npoints dans scheff2 c |
---|
183 | c Declaration real zloc, gam, beami, brhof, radian , c |
---|
184 | c twopi, z0 dans out6dp c |
---|
185 | c Dans les appels a scheff ajouter 2 variables mt c |
---|
186 | c et wtsauv c |
---|
187 | c Modification de swap pour /pcord/ c |
---|
188 | c Suppression des instructions c |
---|
189 | c dimension ac(14) cc(9) bc(1+(1+lmx)*13 c |
---|
190 | c equivalence (ac(1),pi) (bc(1),nel) (cc(1),wt0) c |
---|
191 | c Ajouter mt aux arguments de la subroutine output c |
---|
192 | c Debug version modif main pardyn ouvrir trwave c |
---|
193 | c new common debug , new possibility for ip c |
---|
194 | c negative value means keep file parmdebug c |
---|
195 | c common/debug/idebug et tous les idebug c |
---|
196 | c octobre 2002 augmentation du nombre de cartes cell numc 20-->60 c |
---|
197 | c c |
---|
198 | c 29/10/2002 nom de l'exe sous OSF et HP-UX : parmela c |
---|
199 | c c |
---|
200 | c nom de l'exe sous OSF et HP-UX : parmela c |
---|
201 | c c |
---|
202 | c 19/04/2007 introduction de la polarisation ksi1,ksi2,ksi3 c |
---|
203 | c modification lecture input4 c |
---|
204 | c modification out6dp et sortie c |
---|
205 | c 18 variables dans ndes1 et ndes2 dans out6 c |
---|
206 | c c |
---|
207 | c Last revision 14/11/2007 mise sous Linux c |
---|
208 | c utilisation du debugger c |
---|
209 | c gdb /exp/sera/mouton/Parmela/V503/Linux/parmelalinux c |
---|
210 | c pour cela faire cp toto.in parmin avant de lancer gdb c |
---|
211 | c c |
---|
212 | c 14/11/2007 modification des open de file.f, pour utilisation c |
---|
213 | c du disque local de la machine /tmp c |
---|
214 | c dans le parmela.sh il faut : c |
---|
215 | c mv /tmp/parmdesz toto.desz c |
---|
216 | c Last revision 18/06/2008 c |
---|
217 | c un seul source de parmela v503 compatible OSF LINUX c |
---|
218 | c un seul nom d'exe parmela c |
---|
219 | c c |
---|
220 | c Release octobre 2009 V5.03 c |
---|
221 | c modification du calcul de dzz dans fieldlal c |
---|
222 | c dzz=z(2)-z(1) ----> dzz=(zmax-zmin)/(nchp-1) c |
---|
223 | c augmentation du nombre de points pour le champ magnetique c |
---|
224 | c nspl 1001 ---> 10001 c |
---|
225 | c nouveau parametre nptcb = 10000 dans param_sz.h c |
---|
226 | c creation d'une nouvelle carte nfieldb avec comme parametre c |
---|
227 | c nptcbu nombre de points que l'on veut utiliser pour le c |
---|
228 | c champ magnetique < ou = 10000 c |
---|
229 | c 04/02/2010 retrait de la carte nfieldb c |
---|
230 | c modification de la 1er carte coil c |
---|
231 | c des tests dans field, solnoid si nb pt champ b > nptcb c |
---|
232 | c c |
---|
233 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
---|
234 | c----------------------------------------------------------------------- |
---|
235 | |
---|
236 | c modif F. Touze, G. Le Meur 15/11/2012 |
---|
237 | c---->program parmela |
---|
238 | subroutine parmela(work_path, lonnom) |
---|
239 | c----------------------------------------------------------------------- |
---|
240 | c label n --> n*100 |
---|
241 | c go to 7000 --> 7008 |
---|
242 | c do 8000 --> 8017 |
---|
243 | c format 9000 --> 9025 |
---|
244 | c----------------------------------------------------------------------- |
---|
245 | c |
---|
246 | include 'param_sz.h' |
---|
247 | include 'var_char.h' |
---|
248 | include 'coordcom.h' |
---|
249 | include 'constcom.h' |
---|
250 | include 'debug.h' |
---|
251 | include 'bfieldcom.h' |
---|
252 | include 'flagcom.h' |
---|
253 | include 'misccom.h' |
---|
254 | include 'ncordscom.h' |
---|
255 | include 'outcom.h' |
---|
256 | include 'phizoutcom.h' |
---|
257 | include 'pipecom.h' |
---|
258 | include 'sccom.h' |
---|
259 | include 'syscom.h' |
---|
260 | include 'tstepcom.h' |
---|
261 | include 'ucom.h' |
---|
262 | include 'upawcom.h' |
---|
263 | |
---|
264 | character*132 work_path |
---|
265 | character*256 nom_input |
---|
266 | character*256 nomfic |
---|
267 | character*256 fileelec |
---|
268 | character*256 ficpoiss |
---|
269 | integer lonnom |
---|
270 | |
---|
271 | character*25 kvers |
---|
272 | common/aswap/acor(7,imaa) |
---|
273 | common/back/bmax,byz0,byzs,pzowmin,pzowmax,iback |
---|
274 | common/cathcur/rcur |
---|
275 | common/chars/title |
---|
276 | common/cout7/corout7(6,imaa),bgs,nelal,nbufl,iout7 |
---|
277 | common/ct/var,rp,dzmin,dzmax,wrm,dwbm,xrmsbg,yrmsbg, |
---|
278 | 1xsmax,xpsmax,ysmax,ypsmax,istat,nebm |
---|
279 | common/ddat/ctime,ddat |
---|
280 | common/errors/eraln(5,0:lmx) |
---|
281 | common/image/jj,ij |
---|
282 | common/intype/intype ! mcd |
---|
283 | common/jones/ajones,zjones,zcath(imaa) ! mcd |
---|
284 | common/out6/nbphase,nbremesh |
---|
285 | common/outbuf/outcor(8,imaa*imb) |
---|
286 | common/outpar/noutp,outpar(0:lmx,10),moutp(10) |
---|
287 | common/part/nparticle(imaa) |
---|
288 | common/pnum/ipnum(imaa) |
---|
289 | common/save/zrsav |
---|
290 | common/sbload1/tsbparm(5,50) |
---|
291 | common/tapes/ifl(2) |
---|
292 | common/version/kvers |
---|
293 | common /wtstep/wtstep,nstep ! mcd |
---|
294 | c-------------------------------------------------------------------------- |
---|
295 | |
---|
296 | character Istart_date*9,Iend_date*9,Istart_time*8,Iend_time*8 |
---|
297 | real*4 time_begin |
---|
298 | |
---|
299 | c_with cernlib |
---|
300 | c call timest(1.0e9) |
---|
301 | c call timed(t0) |
---|
302 | c |
---|
303 | |
---|
304 | write(*,*) ' Go Parmela' |
---|
305 | |
---|
306 | |
---|
307 | c_for ibm |
---|
308 | c time_begin = secnds(0.0) |
---|
309 | |
---|
310 | call datime(Istart_date,Istart_time) |
---|
311 | c* |
---|
312 | c---set constants |
---|
313 | c |
---|
314 | kvers='V503 on LINUX' |
---|
315 | c kvers='V503 on OSF' |
---|
316 | c |
---|
317 | pi=4.*atan(1.) |
---|
318 | twopi=2.*pi |
---|
319 | radian=pi/180. |
---|
320 | clight=29979.29 ! 10.e6 cm/s |
---|
321 | erest=.511 ! MeV |
---|
322 | brhof=1704.51 ! gauss.cm moc/e |
---|
323 | c brhof = (9.109389E-31*299792458/1.602177E-19)*(100)*(10000) !1704.5 |
---|
324 | zlimit=1000000. ! cm |
---|
325 | c---set flags |
---|
326 | idebug=0 ! for debuging in subroutine ouvrir |
---|
327 | ifld=0 ! flag for external electric field |
---|
328 | poiflag=.false. ! flag for external magnetic field |
---|
329 | nel=0 ! number of element |
---|
330 | ip=0 ! flag of output type |
---|
331 | isc=0 ! flag for space charge calculation |
---|
332 | istat=0 ! flag for statistics |
---|
333 | ifl(1)=0 |
---|
334 | ifl(2)=0 |
---|
335 | iback=0 |
---|
336 | ltype=1 ! by default disk and washer cells linac |
---|
337 | ifoclal=0 ! bm bzlal |
---|
338 | npipe=0 ! mcd |
---|
339 | nstep=0 ! mcd |
---|
340 | wtstep=0. ! mcd |
---|
341 | ajones=0. ! mcd |
---|
342 | nout=0 ! mcd |
---|
343 | do 8001 n=1,10 ! mcd |
---|
344 | moutp(n)=0 ! mcd |
---|
345 | do 8001 i=1,lmx ! mcd |
---|
346 | outpar(i,n)=0. ! mcd |
---|
347 | 8001 continue ! mcd |
---|
348 | icfield = 0 ! cbm |
---|
349 | ictrw= 0 ! cbm |
---|
350 | nbphase=0 ! cbm nombre de phase en sortie |
---|
351 | nbremesh=0 ! cbm nombre de remesh |
---|
352 | c----set polarization to zero |
---|
353 | do 7999 n=1,imaa |
---|
354 | ksi1(n)=0. |
---|
355 | ksi2(n)=0. |
---|
356 | ksi3(n)=0. |
---|
357 | 7999 continue |
---|
358 | c--- |
---|
359 | c--- |
---|
360 | nom_input = nomfic(work_path(1:lonnom),'parmin') |
---|
361 | print *, ' parmela : fichier input : ', nom_input |
---|
362 | open(unit=nin,file=nom_input,status='old', |
---|
363 | * access='sequential',form='formatted') |
---|
364 | |
---|
365 | |
---|
366 | c$$$ open(unit=nin,file='parmin',status='old', |
---|
367 | c$$$ * access='sequential',form='formatted') |
---|
368 | open(unit=nnout,file='parmdebug',status='unknown', |
---|
369 | * access='sequential',form='formatted') |
---|
370 | call ouvrir(lonnom, work_path) |
---|
371 | c--- |
---|
372 | call setwall |
---|
373 | c---set misalignments to zero. |
---|
374 | do 8002 n=1,100 |
---|
375 | do 8002 i=1,5 |
---|
376 | eraln(i,n)=0. |
---|
377 | 8002 continue |
---|
378 | c---read input data |
---|
379 | c---the read(jj) reads one card with a label and all |
---|
380 | c---subsequent cards without labels. the stop subroutine compares the label |
---|
381 | c---with a list of words and assigns a value to jj corresponding to the |
---|
382 | c---position of the word in the dictionary. the word list and corresponding |
---|
383 | c---values of jj are: |
---|
384 | c---(drift ,1), (solenoid ,2), (quad ,3), (bend ,4), |
---|
385 | c---(buncher ,5), (chopper ,6), (cell ,7), (tank ,8), |
---|
386 | c---(trwave ,9), (coil ,10),(run ,11),(input ,12), |
---|
387 | c---(output ,13),(title ,14),(scheff ,15),(zout ,16), |
---|
388 | c---(adjust ,17),(start ,18),(restart ,19),(continue ,20), |
---|
389 | c---(save ,21),(end ,22),(zlimit ,23),(errors ,24), |
---|
390 | c---(change ,25),(rotate ,26),(sbload ,27),(cfield ,28), |
---|
391 | c---(dpout ,29),(cathode ,30),(design ,31),(pipe ,32), |
---|
392 | c---(foclal ,33),(backb ,34),(wiggler ,35),(alpham ,36), |
---|
393 | c---(stat ,37),(poisson ,38),(sextupole ,39) |
---|
394 | c---the subroutine reads numbers from the card in a free format. |
---|
395 | c---the values are stored in vv, and there are nn of them. |
---|
396 | 99 continue |
---|
397 | cbm call cread(jj) |
---|
398 | call cread |
---|
399 | go to (100,100,100,100,100,100,100,100,100,1000, |
---|
400 | 1 1100,1200,1300,1400,1500,1600,1700,1800,1900,2000, |
---|
401 | 2 2100,2200,2300,2400,2500,100,100,2800,100,100, |
---|
402 | 3 3100,3200,3300,3400,100,100,3700,3800,100),jj |
---|
403 | c---after transport element or linac cell |
---|
404 | c---after drift: vv(n)=l,aper,iout |
---|
405 | c---after solenoid: vv(n)=l,aper,iout,h |
---|
406 | c---after quad: vv(n)=l,aper,iout,hp |
---|
407 | c---after chopper: vv(n)=0,aper,iout,freq,phi,dphi |
---|
408 | c---after cathode: vv(n)=0,aper,iout,radius of curvate of cathode |
---|
409 | c---after wiggler: vv(n)=l,aper,iout,periods,steps,b0,ky/kw |
---|
410 | c---after sextupole: vv(n)=l,aper,iout,Bpole/ape**2 |
---|
411 | 100 continue |
---|
412 | if(nn.gt.11.and.jj.ne.7.and.jj.ne.9.and.jj.ne.27)go to 7000 |
---|
413 | if(nel.ge.lmx)go to 7001 |
---|
414 | if(nn.lt.2)go to 7002 |
---|
415 | nel=nel+1 |
---|
416 | ntype(nel)=jj |
---|
417 | do 8003 n=1,nn |
---|
418 | el(n,nel)=0. ! mcd |
---|
419 | if(n.gt.nn)go to 8003 |
---|
420 | el(n,nel)=vv(n) |
---|
421 | 8003 continue |
---|
422 | zloc(nel)=el(1,nel) |
---|
423 | if(nel.gt.1)zloc(nel)=zloc(nel)+zloc(nel-1) |
---|
424 | if (ip.eq.1 .or. ip.eq.2) then |
---|
425 | iout = vv(3) |
---|
426 | write(nnout,9001) nel,jj,vv(1),zloc(nel),vv(2),iout |
---|
427 | 9001 format ( |
---|
428 | . ' element no. [ne] =',t40,i8/ |
---|
429 | . ' element type [ntype(ne)] =',t40,i8/ |
---|
430 | . ' element length [el(1,ne)] =',t40,1pg12.4,' cm'/ |
---|
431 | . ' z coord at downstream end [zloc(ne)] =',t40,g12.4,' cm'/ |
---|
432 | . ' radial aperture at exit [el(2,ne)] =',t40,g12.4,' cm'/ |
---|
433 | . ' output at end of element [el(3,ne)] =',t40,i8) |
---|
434 | if((jj.eq.3) .or. (jj.eq.39)) write(nnout,9002) vv(4) |
---|
435 | 9002 format( |
---|
436 | . ' magnetic field gradient = ',t40,g12.3,' gauss/cm') |
---|
437 | endif |
---|
438 | if(jj.eq.2) ifld=1 |
---|
439 | if(jj.eq.26) go to 2600 |
---|
440 | if(jj.eq.29) go to 99 |
---|
441 | if(jj.eq.27) go to 2700 |
---|
442 | if(jj.eq.30) rcur=vv(4) |
---|
443 | if(jj.eq.35) go to 3500 |
---|
444 | if(jj.eq.36) go to 3600 |
---|
445 | if(jj.eq.39) go to 99 |
---|
446 | go to (99,99,99,400,500,99,700,800,900), jj |
---|
447 | c---after bend: vv(n)=l,aper,iout,wr,alpha,beta1,beta2,psi1,psi2,r1,r2 |
---|
448 | 400 continue |
---|
449 | do 8004 n=5,9 |
---|
450 | el(n,nel)=el(n,nel)*radian |
---|
451 | 8004 continue |
---|
452 | go to 99 |
---|
453 | c---after buncher: vv(n)=0.,aper,iout,wmax,freq,phi,w |
---|
454 | 500 continue |
---|
455 | w=w0 |
---|
456 | if(nn.ge.7)w=el(7,nel) |
---|
457 | g=w/erest |
---|
458 | bg=sqrt(g*(2.+g)) |
---|
459 | cay=twopi*el(5,nel)/(bg*clight) |
---|
460 | el(8,nel)=cay |
---|
461 | go to 99 |
---|
462 | c---after trwave; |
---|
463 | 900 continue |
---|
464 | ictrw = ictrw+1 |
---|
465 | call trwfield |
---|
466 | go to 99 |
---|
467 | c---after cell: vv(n)=cl,aper,iout,phi,e0,nc,(dwtmax) |
---|
468 | 700 continue |
---|
469 | el(5,nel)=.01*el(5,nel) |
---|
470 | if (nn.lt.8) el(8,nel)=0. |
---|
471 | phideg(nel)=el(4,nel) |
---|
472 | c phi=el(4,nel)*radian this is now defined in cellfld |
---|
473 | c el(9,nel)=sin(phi) and this |
---|
474 | c el(10,nel)=cos(phi) and this also |
---|
475 | c---if unspecified, set maximum step size in cell to 10 degrees. |
---|
476 | if(nn.lt.7)el(7,nel)=10. |
---|
477 | c---get field information for this cell |
---|
478 | cl=el(1,nel) |
---|
479 | if (el(8,nel).ne.0.) cl=2.*cl |
---|
480 | nc=el(6,nel) |
---|
481 | c iel(6,nel)=nc |
---|
482 | if (ip.eq.1 .or. ip.eq.2) then |
---|
483 | write(nnout,9003)vv(4),vv(5),nc,el(7,nel),el(11,nel),el(8,nel) |
---|
484 | 9003 format ( |
---|
485 | . ' phase of cell [phi0] =',t40,g12.3,' degrees'/ |
---|
486 | . ' average axial electric field [el(5,ne)] =',t40,g12.3,' MV/m'/ |
---|
487 | . ' cell no. [nc] =',t40,i8/ |
---|
488 | . ' max step size [dwtmax] =',t40,g12.3,' degrees'/ |
---|
489 | . ' axial static magnetic field =',t40,g12.3,' gauss'/ |
---|
490 | . ' cell symmetry [el(8,ne)] =',t40,g12.3) |
---|
491 | if (el(8,nel).eq.0.) write(nnout,9004) |
---|
492 | 9004 format (' cell is symmetric about midplane in z') |
---|
493 | if (el(8,nel).gt.0.) write(nnout,9005) |
---|
494 | 9005 format (' half cell, with zero field at entrance') |
---|
495 | if (el(8,nel).lt.0.) write(nnout,9006) |
---|
496 | 9006 format (' half cell, with zero field at exit') |
---|
497 | endif |
---|
498 | call cellfld(cl,nc) |
---|
499 | go to 99 |
---|
500 | c---after tank: tl,aper,iout,e0t,nc,p,phi |
---|
501 | 800 continue |
---|
502 | el(8,nel)=el(1,nel)/el(5,nel) |
---|
503 | c---convert e0t from mv/m to mv/cm |
---|
504 | el(4,nel)=.01*el(4,nel) |
---|
505 | c---calculate beta and k**2 for the tank |
---|
506 | beta=2.*el(8,nel)/(el(6,nel)*wavel) |
---|
507 | cl=el(8,nel) |
---|
508 | caysq=((wavel/(2.*cl))**2 - 1.)*(twopi/wavel)**2 |
---|
509 | el(9,nel)=beta |
---|
510 | el(10,nel)=caysq |
---|
511 | go to 99 |
---|
512 | c---after field: vv parameters depend on field subroutine |
---|
513 | 1000 continue |
---|
514 | icfield=icfield+1 |
---|
515 | call field |
---|
516 | c---define field structure after coil |
---|
517 | c--- vv(n)=z,r,ni,(zbeg,zend,nptcbu,order) |
---|
518 | c---in ( ) required on first coil card ignored on rest |
---|
519 | go to 99 |
---|
520 | c---after errors. generate errors in geners |
---|
521 | 2400 continue |
---|
522 | call geners |
---|
523 | go to 99 |
---|
524 | c---after run: vv(n)=irun,ip,freq,z0,w0,ltype |
---|
525 | 1100 continue |
---|
526 | c process run card ! mcd |
---|
527 | if (nn.ge.1) irun=vv(1) ! mcd |
---|
528 | c--debug |
---|
529 | c-- Si ip negatif sur carte run on a la version de debug |
---|
530 | if (nn.ge.2) then |
---|
531 | if(vv(2).lt.0.) then |
---|
532 | idebug=1 |
---|
533 | write(*,*) '-----------------------------------' |
---|
534 | write(*,*) 'You have ask the version with debug' |
---|
535 | write(*,*) '-----------------------------------' |
---|
536 | write(nnout,*) '-----------------------------------' |
---|
537 | write(nnout,*) 'You have ask the version with debug' |
---|
538 | write(nnout,*) '-----------------------------------' |
---|
539 | else |
---|
540 | endif |
---|
541 | ip=abs(vv(2)) |
---|
542 | endif |
---|
543 | c--debug |
---|
544 | if (nn.ge.3) freq=vv(3) ! mcd |
---|
545 | if (nn.lt.4) go to 99 |
---|
546 | if (nn.ge.4) z0=vv(4) ! mcd |
---|
547 | if (nn.lt.5) go to 99 |
---|
548 | if (nn.ge.5) w0=vv(5) ! mcd |
---|
549 | if (nn.ge.6) ltype=vv(6) ! mcd |
---|
550 | if(freq.gt.0.) wavel=clight/freq |
---|
551 | write(nnout,9007)irun |
---|
552 | 9007 format (/' '/ |
---|
553 | . ' run no. =',t40,i8/) |
---|
554 | if((ip.gt.4).and.(ip.ne.999)) then |
---|
555 | ip=1 |
---|
556 | write(nnout,9008)ip |
---|
557 | 9008 format( |
---|
558 | .' print option modified =',t40,i8) |
---|
559 | else |
---|
560 | write(nnout,9009)ip |
---|
561 | 9009 format( |
---|
562 | . ' print option =',t40,i8) |
---|
563 | endif |
---|
564 | c--- |
---|
565 | if (ip.eq.3) then |
---|
566 | open(unit=nuch,file='parmtrch',access='sequential', |
---|
567 | * status='unknown') |
---|
568 | write(nuch,*)' nc sym r z ez er bt ' |
---|
569 | endif |
---|
570 | if (ip.eq.999) open(unit=nimp,file='parmimp',status='unknown') |
---|
571 | c--- |
---|
572 | write(nnout,9010)freq,wavel,z0,w0,ltype |
---|
573 | 9010 format( |
---|
574 | . ' rf frequency =',t40,f9.2,' MHz'/ |
---|
575 | . ' wavelength =',t40,1pg12.3,' cm'/ |
---|
576 | , ' z0 =',t40,g12.3,' cm'/ |
---|
577 | , ' reference particle energy =',t40,g12.3,' MeV'/ |
---|
578 | . ' linac type =',t40,i8) |
---|
579 | do 8005 i=1,4 |
---|
580 | corout7(i,1)=0. |
---|
581 | cor(i,1)=0. |
---|
582 | 8005 continue |
---|
583 | cor(5,1)=z0 |
---|
584 | corout7(5,1)=z0 |
---|
585 | corout7(6,1)=w0 |
---|
586 | g=w0/erest |
---|
587 | gam(1)=g+1 |
---|
588 | cor(6,1)=sqrt(g*(2.+g)) |
---|
589 | cord(6,1)=cor(6,1) |
---|
590 | if (npoints.lt.1)npoints=1 |
---|
591 | ngood=npoints |
---|
592 | nparticle(1)=1 |
---|
593 | phizero(1)=0. |
---|
594 | if(nn.lt.6)go to 99 |
---|
595 | go to 99 |
---|
596 | c---after input: vv parameters depend on input subroutine |
---|
597 | 1200 continue |
---|
598 | call input(lonnom, work_path) |
---|
599 | ngood=npoints |
---|
600 | do 8006 i=1,ngood |
---|
601 | nparticle(i)=i |
---|
602 | ipnum(i)=i |
---|
603 | 8006 continue |
---|
604 | go to 99 |
---|
605 | c---after output: vv parameters stored in optcon array |
---|
606 | c---after output: vv parameters stored in outpar array |
---|
607 | c |
---|
608 | 1300 continue |
---|
609 | noutp=vv(1)+1 |
---|
610 | if(noutp.gt.10 .or. noutp.lt.1) then |
---|
611 | write(ndiag,9011) |
---|
612 | 9011 format(' invalid OUTPUT type') |
---|
613 | go to 99 |
---|
614 | endif |
---|
615 | moutp(noutp)=nn |
---|
616 | do 8007 n=1,nn |
---|
617 | optcon(n)=vv(n) |
---|
618 | outpar(n,noutp)=vv(n) |
---|
619 | 8007 continue |
---|
620 | if(vv(1).gt.1..and.vv(1).lt.5.)then |
---|
621 | imbf=imb-1 |
---|
622 | else |
---|
623 | imbf=imb |
---|
624 | endif |
---|
625 | go to 99 |
---|
626 | c---after title: read title information from next card |
---|
627 | 1400 continue |
---|
628 | read(nin,9012)title |
---|
629 | 9012 format(a80) |
---|
630 | write(nnout,*)title |
---|
631 | go to 99 |
---|
632 | c---after scheff: vv(1)=beam current; vv(2)= program number |
---|
633 | c--- rest depend on space-charge subroutines |
---|
634 | 1500 continue |
---|
635 | do 8008 n=1,12 |
---|
636 | sce(n)=0. |
---|
637 | if (n.le.nn) sce(n)=vv(n) |
---|
638 | 8008 continue |
---|
639 | beami=vv(1) |
---|
640 | iprog=vv(2) |
---|
641 | if(isc.eq.0) then |
---|
642 | if(iprog.eq.1) call scheff1(0.,-99,0.) |
---|
643 | if(iprog.eq.2) call scheff2(0.,-99,0.) |
---|
644 | if(iprog.eq.3) call scheff3(0.,-99,0.) |
---|
645 | isc=1 |
---|
646 | else ! isc |
---|
647 | endif ! isc |
---|
648 | go to 99 |
---|
649 | c---after zout. print locations of elements |
---|
650 | 1600 continue |
---|
651 | if(ifoclal.eq.1) then |
---|
652 | call zoutlal |
---|
653 | else |
---|
654 | call zout |
---|
655 | endif |
---|
656 | go to 99 |
---|
657 | c---after adjust: call arbitrary subroutine named adjust |
---|
658 | 1700 continue |
---|
659 | call adjust |
---|
660 | go to 99 |
---|
661 | c---after change. call change subroutine |
---|
662 | 2500 continue |
---|
663 | call change |
---|
664 | go to 99 |
---|
665 | c---after rotate |
---|
666 | 2600 continue |
---|
667 | el(4,nel)=el(4,nel)*radian |
---|
668 | go to 99 |
---|
669 | c--- after sbload paramaters are l,aperature,iout,zlen,ncells,a1,a2,a3 |
---|
670 | c--- ,a4,a5,ta1,ta2,ta3,ta4,ta5. if ncells is negative then the paramaters |
---|
671 | c--- ta1 to ta4 are required. ta5 is optional. if ncells is negative the |
---|
672 | c--- transvers wake will be applied also. |
---|
673 | 2700 continue |
---|
674 | if(nn.gt.11)then |
---|
675 | numsb=numsb+1 |
---|
676 | tsbparm(1,numsb)=vv(11) |
---|
677 | el(11,nel)=numsb |
---|
678 | tsbparm(2,numsb)=vv(12) |
---|
679 | tsbparm(3,numsb)=vv(13) |
---|
680 | tsbparm(4,numsb)=vv(14) |
---|
681 | if(el(5,nel).lt.0..and.nn.lt.14)go to 7002 |
---|
682 | if(nn.gt.14)tsbparm(5,numsb)=vv(15) |
---|
683 | endif |
---|
684 | go to 99 |
---|
685 | c---after start: start dynamics calculations |
---|
686 | 1800 continue |
---|
687 | if(nn.lt.5)go to 7002 |
---|
688 | wt0=vv(1) |
---|
689 | c--- For input type 10 and 11, WTSTEP is the half width of the laser pulse. |
---|
690 | c The clock must be set back by this amount to generate the time |
---|
691 | c distribution of photoelectrons. |
---|
692 | wt0=wt0-wtstep |
---|
693 | dwt=vv(2) |
---|
694 | nsteps=vv(3) |
---|
695 | nsc=vv(4) |
---|
696 | nout=vv(5) |
---|
697 | go to 7003 |
---|
698 | c---after restart. set the starting time (phase), initialize cord |
---|
699 | c---from cor, and find largest element number (len) that any particle is in. |
---|
700 | 1900 if(nn.lt.4)go to 7002 |
---|
701 | dwt=vv(1) |
---|
702 | nsteps=vv(2) |
---|
703 | nsc=vv(3) |
---|
704 | nout=vv(4) |
---|
705 | c ----------------- |
---|
706 | if(jj.eq.19) then |
---|
707 | npoints=-vv(6) |
---|
708 | endif |
---|
709 | c ------------------------- |
---|
710 | if(nn.gt.5.and.(vv(6).gt.0..and.vv(6).le.26))then |
---|
711 | i=vv(6)-1 |
---|
712 | call orest1(i) |
---|
713 | else |
---|
714 | call orest2 |
---|
715 | endif |
---|
716 | read(nsav) cor,wt0,ngs |
---|
717 | close(nsav) |
---|
718 | ngood=ngs |
---|
719 | C ---------------------------- |
---|
720 | IF(NGS.NE.NPOINTS) then |
---|
721 | call appendparm |
---|
722 | print *,'ngs=',ngs,'npoints=',npoints |
---|
723 | stop ' Abnormal stop ngs.ne.npoints main ' |
---|
724 | endif |
---|
725 | C ------------------------------- |
---|
726 | 7003 continue |
---|
727 | if (ip.eq.1 .or. ip.eq.2) |
---|
728 | . write (nnout,9013) wt0,dwt,nsteps,vv(5),nsc,nout |
---|
729 | 9013 format ( |
---|
730 | . ' starting phase =',t40,1pg12.3,' degrees'/ |
---|
731 | . ' integration step size =',t40,g12.3,' degrees'/ |
---|
732 | . ' no. of steps in integration =',t40,i8/ |
---|
733 | . ' shift of reference particle =',t40,g12.3,' cm'/ |
---|
734 | . ' call SCHEFF every',t40,i8,' step(s)'/ |
---|
735 | . ' call OUTPUT every',t40,i8,' step(s)'/) |
---|
736 | wt=wt0 |
---|
737 | do 8009 i=1,maxbuf |
---|
738 | neb(i)=0 |
---|
739 | nbuf(i)=0 |
---|
740 | 8009 continue |
---|
741 | c---if nn.gt.4, the 5th parameter is dist to shift ref particle |
---|
742 | c only for restart. not for start |
---|
743 | if (jj.eq.19.and.nn.gt.4) then |
---|
744 | cor(5,1)=cor(5,1)+vv(5) |
---|
745 | endif |
---|
746 | c -------------------------------------------------- |
---|
747 | c--------- vv(5)= dz |
---|
748 | if(jj.eq.19.and.vv(5).ne.0.) then |
---|
749 | cor(1,1)=0. ! for energy change of ref part |
---|
750 | cor(2,1)=0. |
---|
751 | cor(3,1)=0. |
---|
752 | cor(4,1)=0. |
---|
753 | print *,'vv(7) =',vv(7) |
---|
754 | cor(6,1)=sqrt((1.+vv(7)/0.511)**2-1.) |
---|
755 | endif |
---|
756 | c --------------------------------------------------- |
---|
757 | do 8010 n=1,ngood |
---|
758 | do 8011 i=1,6 |
---|
759 | cord(i,n)=cor(i,n) |
---|
760 | 8011 continue |
---|
761 | gam(n)=sqrt(1.+cord(2,n)**2+cord(4,n)**2+cord(6,n)**2) |
---|
762 | c---find which element each particle is in |
---|
763 | ne=0 |
---|
764 | z=cord(5,n) |
---|
765 | if(z.lt.0.)go to 7004 |
---|
766 | do 8012 i=1,nel |
---|
767 | ne=i |
---|
768 | if(z.lt.zloc(i))go to 7004 |
---|
769 | 8012 continue |
---|
770 | ne=nel+1 |
---|
771 | 7004 cord(7,n)=ne |
---|
772 | c if generating from a photocathode, put ne=0 |
---|
773 | if(jj.eq.19) go to 8010 |
---|
774 | if (intype.eq.10 .or. intype.eq.11 .or. intype.eq.15) |
---|
775 | . cord(7,n)=0. |
---|
776 | 8010 continue |
---|
777 | c---if nn.gt.4, the 5th parameter is dist to shift ref particle |
---|
778 | c only for restart. not for start |
---|
779 | c---start dynamics calculation |
---|
780 | call pardyn |
---|
781 | go to 99 |
---|
782 | c---after continue: used to change step size and continue computations |
---|
783 | 2000 continue |
---|
784 | if(nn.lt.4)go to 7002 |
---|
785 | dwt=vv(1) |
---|
786 | nsteps=vv(2) |
---|
787 | nsc=vv(3) |
---|
788 | nout=vv(4) |
---|
789 | c---if nn.gt.4, the 5th parameter is dist to shift ref particle |
---|
790 | if (nn.gt.4)cord(5,1)=cord(5,1)+vv(5) |
---|
791 | call pardyn |
---|
792 | go to 99 |
---|
793 | c---after save |
---|
794 | c---save cord in cor |
---|
795 | 2100 continue |
---|
796 | ngs=ngood |
---|
797 | do 8013 n=1,ngood |
---|
798 | do 8013 i=1,6 |
---|
799 | cor(i,n)=cord(i,n) |
---|
800 | 8013 continue |
---|
801 | write(nnout,*) 'zr of save ',zrsav |
---|
802 | wt0=wt |
---|
803 | if(nn.ge.1.and.(vv(1).gt.0..and.vv(1).le.26))then |
---|
804 | i=vv(1)-1 |
---|
805 | call osave1(i) |
---|
806 | if(ifl(2).ne.0) then |
---|
807 | close(nsav) |
---|
808 | ifl(2)=0 |
---|
809 | endif |
---|
810 | if(idebug.eq.0) then |
---|
811 | close(unit=nnout,status='delete') |
---|
812 | else |
---|
813 | close(unit=nnout,status='keep') |
---|
814 | endif |
---|
815 | close(nnout) |
---|
816 | open(unit=nnout,file='parmout',access='sequential', |
---|
817 | * status='unknown',form='formatted') |
---|
818 | 7005 continue |
---|
819 | read(nnout,'(10a8)',err=7006)tem |
---|
820 | go to 7005 |
---|
821 | 7006 continue |
---|
822 | backspace (nnout) |
---|
823 | else |
---|
824 | open (unit=nsav,file='savecor',access='sequential', |
---|
825 | * form='unformatted',status='new') |
---|
826 | endif |
---|
827 | write(nsav) cor,wt0,ngs,bgs |
---|
828 | close(nsav) |
---|
829 | go to 99 |
---|
830 | c---after zlimit: lower limit on lagging particles |
---|
831 | 2300 continue |
---|
832 | zlimit=vv(1) |
---|
833 | go to 99 |
---|
834 | c---after cfield: read fields from sfo7 as written on tape7 |
---|
835 | c--- first parameter on cfield is starting cell number |
---|
836 | c--- next parameter is ending cell number and third parameter is |
---|
837 | c--- increment of cell number as in a do loop |
---|
838 | c---line after cfield is filename of sfo7 data |
---|
839 | 2800 continue |
---|
840 | C modif glm 29/10/2012 read(nin,'(a12)')filesuper |
---|
841 | read(nin,*) filesuper |
---|
842 | print *, ' j''ai lu le fichier champ : ', filesuper |
---|
843 | write(nnout,'(1x,a12)')filesuper |
---|
844 | if(nn.lt.1)then |
---|
845 | nc=1 |
---|
846 | ne=1 |
---|
847 | ns=1 |
---|
848 | elseif(nn.eq.1)then |
---|
849 | nc=vv(1) |
---|
850 | ne=nc |
---|
851 | ns=1 |
---|
852 | elseif(nn.eq.2)then |
---|
853 | nc=vv(1) |
---|
854 | ne=vv(2) |
---|
855 | ns=1 |
---|
856 | elseif(nn.ge.3)then |
---|
857 | nc=vv(1) |
---|
858 | ne=vv(2) |
---|
859 | ns=vv(3) |
---|
860 | endif |
---|
861 | |
---|
862 | c$$$ call cellfeld(filesuper,nc,ne,ns) |
---|
863 | fileelec = nomfic(work_path(1:lonnom),filesuper) |
---|
864 | print *, ' j''envoie lu le fichier champ : ', fileelec |
---|
865 | call cellfeld(fileelec,nc,ne,ns) ! modif glm 16/11/2012 |
---|
866 | go to 99 |
---|
867 | c---after design |
---|
868 | 3100 continue |
---|
869 | call design |
---|
870 | go to 99 |
---|
871 | c---after pipe, vv(n) = radius,zlow,zhigh |
---|
872 | 3200 continue |
---|
873 | if (npipe.eq.100) then |
---|
874 | write(nnout,9014) |
---|
875 | 9014 format (' Too many PIPE cards') |
---|
876 | go to 99 |
---|
877 | endif |
---|
878 | if (vv(1).le.0.) then |
---|
879 | write(nnout,9015) |
---|
880 | 9015 format (' Bad radius on PIPE card') |
---|
881 | go to 99 |
---|
882 | endif |
---|
883 | npipe = npipe + 1 |
---|
884 | rpipe(npipe) = vv(1) |
---|
885 | zlpipe(npipe) = vv(2) |
---|
886 | zhpipe(npipe) = vv(3) |
---|
887 | if (ip.eq.1 .or. ip.eq.2) |
---|
888 | .write(nnout,9016) npipe,vv(1),vv(2),vv(3) |
---|
889 | 9016 format ( |
---|
890 | . ' Beam pipe for point-by-point image calculation'/ |
---|
891 | . ' pipe no. = ',t40,i8/ |
---|
892 | . ' pipe radius = ',t40,1pg12.3,' cm'/ |
---|
893 | . ' zlow = ',t40,g12.3,' cm'/ |
---|
894 | . ' zhigh = ',t40,g12.3,' cm') |
---|
895 | go to 99 |
---|
896 | c---after foclal vv(n)=zmin,zmax,dzz,npchb,facbz,opt |
---|
897 | 3300 continue |
---|
898 | iopt=vv(6) |
---|
899 | if(iopt.eq.0) then |
---|
900 | ifld=0 |
---|
901 | ifoclal=0 |
---|
902 | read(nin,'(a12)') filebz |
---|
903 | write(nnout,*) 'foclal not used' |
---|
904 | else |
---|
905 | ifld=1 |
---|
906 | ifoclal=1 |
---|
907 | read(nin,'(a12)') filebz |
---|
908 | call fieldlal(filebz) |
---|
909 | endif |
---|
910 | go to 99 |
---|
911 | c---after backb back bombardment |
---|
912 | 3400 continue |
---|
913 | iback=1 |
---|
914 | bmax=vv(1) |
---|
915 | byz0=vv(2) |
---|
916 | byzs=vv(3) |
---|
917 | iope=0 |
---|
918 | if(nn.gt.3) then |
---|
919 | pzowmin=vv(4) |
---|
920 | pzowmax=vv(5) |
---|
921 | else |
---|
922 | pzowmin=0. |
---|
923 | pzowmax=0. |
---|
924 | endif |
---|
925 | if(pzowmax.ne.0.) iope=1 |
---|
926 | call opeback(iope) |
---|
927 | go to 99 |
---|
928 | c---after wiggler l,aper,iout,periods,steps,b0,kx,ky,kw,iprint |
---|
929 | 3500 continue |
---|
930 | if(nn.lt.7) go to 7002 |
---|
931 | xlam=el(1,nel)/el(4,nel) |
---|
932 | cayw=twopi/xlam |
---|
933 | cayy=el(7,nel)*cayw |
---|
934 | if(el(7,nel).le.0. .or. el(7,nel).gt.1.)then |
---|
935 | write(nnout,9018) |
---|
936 | 9018 format(' error: ky/kw on wiggler must be between 0 and 1') |
---|
937 | jj=100 ! for cread |
---|
938 | go to 99 |
---|
939 | endif |
---|
940 | cayx=0. |
---|
941 | cayxsq=cayw**2-cayy**2 |
---|
942 | if(cayxsq.gt.0.) cayx=sqrt(cayxsq) |
---|
943 | el(7,nel)=cayx |
---|
944 | el(8,nel)=cayy |
---|
945 | el(9,nel)=cayw |
---|
946 | if(ip.eq.1 .or. ip.eq.2) |
---|
947 | .write(nnout,9019)el(4,nel),el(5,nel),el(6,nel),el(7,nel) |
---|
948 | 9019 format( |
---|
949 | .' number of wiggler periods =',t40,g12.3/ |
---|
950 | .' step of wiggler =',t40,g12.3/ |
---|
951 | .' wiggler field =',t40,g12.3,' gauss'/ |
---|
952 | .' ratio of Ky to Kw =',t40,g12.3) |
---|
953 | go to 99 |
---|
954 | c---after alpha magnet |
---|
955 | 3600 continue |
---|
956 | if(ip.eq.1 .or. ip.eq.2) |
---|
957 | .write(nnout,9020) el(1,nel),el(4,nel),el(5,nel),el(7,nel) |
---|
958 | 9020 format( |
---|
959 | .' distance travelled by the r.p. =',t40,g12.3,' cm'/ |
---|
960 | .' reference particle energy =',t40,g12.3,' MeV'/ |
---|
961 | .' alpha magnet angle =',t40,g12.3,' Degrees'/ |
---|
962 | .' absolute field gradient =',t40,g12.3,' Gauss/cm') |
---|
963 | go to 99 |
---|
964 | c---after stat |
---|
965 | 3700 continue |
---|
966 | var=vv(1) |
---|
967 | istat=1 |
---|
968 | open(unit=nstat,file='parmstat',form='unformatted', |
---|
969 | * status='unknown') |
---|
970 | write(nstat)ij |
---|
971 | go to 99 |
---|
972 | c---after Poisson |
---|
973 | 3800 continue |
---|
974 | C modif glm 29/10/2012 read (nin,'(a12)')filepois |
---|
975 | read (nin,*)filepois |
---|
976 | write(nnout,'(1x,a12)')filepois |
---|
977 | if(nn.ge.1)then |
---|
978 | zoffset=vv(1) |
---|
979 | else |
---|
980 | zoffset=0. |
---|
981 | endif |
---|
982 | if(nn.ge.2)then |
---|
983 | rmult=vv(2) |
---|
984 | else |
---|
985 | rmult=1. |
---|
986 | endif |
---|
987 | ficpoiss = nomfic(work_path(1:lonnom),filepois) |
---|
988 | call poisson(ficpoiss,zoffset,rmult) |
---|
989 | go to 99 |
---|
990 | c |
---|
991 | cbm retrait de nfieldb 04/02/2010 apres modif de coil |
---|
992 | c--- after nfieldb number of points used for the magnetic field |
---|
993 | c < 10000 |
---|
994 | c 4000 continue |
---|
995 | c nptcbu=vv(1) |
---|
996 | c write(nnout,*) |
---|
997 | c *' Nombre de points maximal pour le champ magnetique ',nptcb |
---|
998 | c write(nnout,*) |
---|
999 | c *' Number of used values for magnetic field',nptcbu |
---|
1000 | c if(nptcbu.gt.nptcb) then |
---|
1001 | c write(nnout,*) |
---|
1002 | c * ' Number of points for magnetic field gt ',nptcb |
---|
1003 | c endif |
---|
1004 | c go to 99 |
---|
1005 | c---after-end output any used buffer |
---|
1006 | 2200 continue |
---|
1007 | c---check output buffers. |
---|
1008 | c---element neo? |
---|
1009 | inb=imaa*imbf/npoints |
---|
1010 | if(inb.gt.maxbuf)inb=maxbuf |
---|
1011 | do 8014 i=1,inb |
---|
1012 | if(nbuf(i).eq.0)go to 8014 |
---|
1013 | neo=neb(i) |
---|
1014 | c---yes. empty buffer. |
---|
1015 | call output(neo,nbuf(i),outcor(1,(i-1)*npoints+1),npoints) |
---|
1016 | 8014 continue |
---|
1017 | do 8015 j=0,nel |
---|
1018 | sum=0. |
---|
1019 | sumi=0. |
---|
1020 | do 8016 i=ngood+1,npoints |
---|
1021 | if(cord(7,i).eq.j)then |
---|
1022 | sum=sum+1 |
---|
1023 | sumi=sumi+(gam(i)-1.)*erest |
---|
1024 | endif |
---|
1025 | 8016 continue |
---|
1026 | if(sum.ne.0.)then |
---|
1027 | sumi=sumi/sum |
---|
1028 | write(nnout,9021) j,sum/npoints,sumi |
---|
1029 | 9021 format(' fraction lost in element',i3,' = ',f7.5, |
---|
1030 | *' average energy =',f9.5,' MeV') |
---|
1031 | endif |
---|
1032 | 8015 continue |
---|
1033 | if(ngood.gt.2)then |
---|
1034 | sumi=0. |
---|
1035 | do 8017 i=1,ngood |
---|
1036 | sumi=sumi+(gam(i)-1.)*erest |
---|
1037 | 8017 continue |
---|
1038 | sumi=sumi/ngood |
---|
1039 | write(nnout,9022) float(ngood)/npoints,sumi |
---|
1040 | 9022 format(' fraction good at end of calculation = ',f7.5, |
---|
1041 | * ' average energy = ',f9.5,' MeV ') |
---|
1042 | endif |
---|
1043 | c----- |
---|
1044 | cbm 18/07/2k if(lolun1.ne.0) call closehisto |
---|
1045 | c----- |
---|
1046 | write(ndiag,*) ' ' |
---|
1047 | write(ndiag,*) ' Statistique sur ce run' |
---|
1048 | write(ndiag,*) ' ' |
---|
1049 | |
---|
1050 | call datime(Iend_date,Iend_time) |
---|
1051 | |
---|
1052 | c-with cernlib |
---|
1053 | c call timed(t1) |
---|
1054 | c CPU_time=t1 |
---|
1055 | c_vms/unix |
---|
1056 | c_vms |
---|
1057 | c write(ndiag,698) |
---|
1058 | c * Istart_date(1:4),Istart_date(5:6),Istart_date(7:8),Istart_time, |
---|
1059 | c * Iend_date(1:4),Iend_date(5:6),Iend_date(7:8),Iend_time, |
---|
1060 | c * secnds(time_begin),CPU_time |
---|
1061 | c 698 format(//,' -- on ALPHA VMS -- LAL --',28('-')/ |
---|
1062 | c * ' Start Execution : ',a4,'/',a2,'/',a2,' ',a8,/ |
---|
1063 | c * ' End Execution : ',a4,'/',a2,'/',a2,' ',a8,/ |
---|
1064 | c * ' Time used for calculation = ',f10.2,' [sec]',/ |
---|
1065 | c * ' CPU Time = ',f10.2,' [sec]',/ |
---|
1066 | c * ' -------------------------------',20('-')) |
---|
1067 | c_unix |
---|
1068 | if (idebug.eq.0) then |
---|
1069 | else |
---|
1070 | write(ndiag,*) ' ---- Version de DEBUG ----- ' |
---|
1071 | endif |
---|
1072 | write(ndiag,698) kvers, |
---|
1073 | * Istart_date,Istart_time, |
---|
1074 | * Iend_date, Iend_time |
---|
1075 | c * ,secnds(time_begin),CPU_time |
---|
1076 | 698 format(//,' -- LAL -- ',a,16('-')/ |
---|
1077 | * ' Start Execution : ',a9,' ',a8,/ |
---|
1078 | * ' End Execution : ',a9,' ',a8,/ |
---|
1079 | * ' -------------------------------',20('-')) |
---|
1080 | c les 2 lignes supprimees apres suppression de la cernlib |
---|
1081 | c retrait des appel a la fonction secnds et a la subroutine timed |
---|
1082 | c * ' Time used for calculation = ',f10.2,' [sec]',/ |
---|
1083 | c * ' CPU Time = ',f10.2,' [sec]',/ |
---|
1084 | |
---|
1085 | c_vms/unix |
---|
1086 | c------------------------------------------------------------------------ |
---|
1087 | call appendparm |
---|
1088 | if(idebug.eq.0) then |
---|
1089 | close(unit=nnout,status='delete') |
---|
1090 | else |
---|
1091 | close(unit=nnout,status='keep') |
---|
1092 | write(*,*) ' ' |
---|
1093 | write(*,*) ' CLOSE statement keep file parmdebug ' |
---|
1094 | write(*,*) ' ' |
---|
1095 | endif |
---|
1096 | print *, ' parmela normal end' |
---|
1097 | stop ' normal end' |
---|
1098 | c---error messages |
---|
1099 | 7000 continue |
---|
1100 | write(nnout,9023) |
---|
1101 | 9023 format(' too many parameters on element card') |
---|
1102 | jj=100 |
---|
1103 | go to 99 |
---|
1104 | 7001 continue |
---|
1105 | write(nnout,9024) lmx |
---|
1106 | 9024 format(' too many elements max is ',i4) |
---|
1107 | go to 2200 |
---|
1108 | 7002 write(nnout,9025) |
---|
1109 | 9025 format(' not enough parameters') |
---|
1110 | jj=100 |
---|
1111 | go to 99 |
---|
1112 | end |
---|
1113 | |
---|
1114 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
---|
1115 | subroutine cellfld(cl,nc) |
---|
1116 | c cell initialization |
---|
1117 | c---------------------------------------------------------------------- |
---|
1118 | save |
---|
1119 | c |
---|
1120 | include 'param_sz.h' |
---|
1121 | include 'cfldscom.h' |
---|
1122 | include 'constcom.h' |
---|
1123 | include 'misccom.h' |
---|
1124 | include 'syscom.h' |
---|
1125 | include 'ucom.h' |
---|
1126 | c |
---|
1127 | common/jones/ajones,zjones,zcath(imaa) |
---|
1128 | common/fcoef/ncoef(100),noff(100),ncer(100) |
---|
1129 | c-------------------------------------------------------------------------- |
---|
1130 | c* |
---|
1131 | hcll(nc)=.5*cl |
---|
1132 | beta=2.*cl/wavel |
---|
1133 | cellfreq(nc)=freq |
---|
1134 | lctype=ltype |
---|
1135 | if(nn.ge.9)then |
---|
1136 | if(vv(9).ne.0.)cellfreq(nc)=vv(9) |
---|
1137 | beta=beta*cellfreq(nc)/freq |
---|
1138 | endif |
---|
1139 | if(nn.ge.10.and.vv(10).ne.0.)lctype=vv(10) |
---|
1140 | phi=el(4,nel)*radian |
---|
1141 | el(9,nel)=sin(phi) |
---|
1142 | el(10,nel)=cos(phi) |
---|
1143 | el(4,nel)=cellfreq(nc) |
---|
1144 | do 30 i = 1,numcf |
---|
1145 | fc(i,nc) = 0. |
---|
1146 | 30 continue |
---|
1147 | if (ip.eq.1 .or. ip.eq.2) then |
---|
1148 | write(nnout,20) lctype,cellfreq(nc),hcll(nc) |
---|
1149 | 20 format ( |
---|
1150 | . ' cell type =',t40,i8/ |
---|
1151 | . ' cell frequency =',t40,1pg12.3,' MHz'/ |
---|
1152 | . ' "half-cell" length =',t40,g12.3,' cm') |
---|
1153 | if (nel.eq.1.and.ajones.gt.0.) then |
---|
1154 | write(nnout,21) |
---|
1155 | 21 format (' rf fields from analytic form of M. Jones') |
---|
1156 | return |
---|
1157 | endif |
---|
1158 | endif |
---|
1159 | if(nn.gt.11) then |
---|
1160 | np=nn |
---|
1161 | ncoef(nel)=vv(12) |
---|
1162 | noff(nel)=vv(13) |
---|
1163 | ncer(nel)=vv(14) |
---|
1164 | nmax=14+ncoef(nel) |
---|
1165 | if(np.gt.nmax)np=nmax |
---|
1166 | do 10 i=15,nmax |
---|
1167 | fc(i-14,nc)=vv(i) |
---|
1168 | 10 continue |
---|
1169 | if (ip.eq.1 .or. ip.eq.2) then |
---|
1170 | write(nnout,39) ncoef(nel),noff(nel) |
---|
1171 | 39 format (' No. of Fourier coefs =',t40,i8/ |
---|
1172 | . ' odd/even cosines =',t40,i8) |
---|
1173 | write(nnout,40) nmax |
---|
1174 | 40 format (' Using vv(15 to ',i2,') as Fourier coefficents of' |
---|
1175 | . ,' this cell') |
---|
1176 | write(nnout,41) ncer(nel) |
---|
1177 | 41 format (' No. of Fourier coefs used in er calculation =', |
---|
1178 | . t50,i8) |
---|
1179 | endif |
---|
1180 | go to 100 |
---|
1181 | endif |
---|
1182 | c get fourier coefficients for special cells from tables |
---|
1183 | ncoef(nel) = 14 |
---|
1184 | noff(nel) = 1 |
---|
1185 | ncer(nel) = 14 |
---|
1186 | iflag(nc)=1 |
---|
1187 | if(iflag(nc).ne.0)return |
---|
1188 | if(lctype.eq.8.and.abs(hcll(nc)-2.21).lt..01)beta=1.02 |
---|
1189 | if(lctype.eq.6)call rtmfc(beta,fc(1,nc)) |
---|
1190 | if(lctype.eq.2)call sccfc(beta,fc(1,nc)) |
---|
1191 | if(lctype.ne.1.and.lctype.ne.2)call genfc(beta,fc(1,nc)) |
---|
1192 | return |
---|
1193 | c make lookup table for rf field |
---|
1194 | 100 continue |
---|
1195 | call genfld(nc) |
---|
1196 | return |
---|
1197 | end |
---|
1198 | c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* |
---|