Brian, This is a simplified version of the program to calculate the trajectories of atomic beams inside of a sextupole magnet. This version has a fixed value for the magnetic moment of the nucleus. The second version of the program which I will also send (ATOMSTATE.FOR) allows the nuclear magnetic moment to be redefined. The code is said to be self documented and includes a 'built-in' dataset - which you will have to customize. The student responsible for this code is Ron Sayer - a UNC- Chapel Hill student. His Bitnet address here at TUNL is "SAYER@TUNL". Feel free to contact him with questions. - Chris Westerfeldt CWEST@TUNL.BITNET ----------------------------cut here-------------------------------------- C C C A T O M S O R C C C C ************** INTRODUCTION **************** C C This is a numerical simulation of particle trajectories in c an atomic beam source, specifically written for the atomic beam c section of the Polarized Ion Source at Triangle Universities c Nuclear Laboratory, but should be fairly adaptable to other set-ups. C Hydrogen is ejected from a nozzle cooled to low temperatures, c estimated to be about 32K. This temperature determines the velocity c distribution of the outgoing particles. The program uses analytic c solutions in its calculations of particle trajectories, each of which c begins with a different axial velocity(VZ), y-position(YI), c y-angle(ALPY), and x-angle(ALPX). Due to the cylindrical symmetry c of the system, x-position and negative x-angles are not needed. c The code can handle any arbitrary velocity distribution for VZ, c inserted in the form of velocity weighting parameters. The peak c of the distribution is called PROBV. The program loops through nine c different velocities in increments of 10000 cm/s. Angular and radial c weightings can likewise be inserted, however none of these distributions c can be measured easily. c This program assumes a nozzle, collimator (called a skimmer), c sextupole 1 (tapered), sextupole 2 (untapered), another collimator (the c iron plate at the entrance to the ionizer region), and the ionizing c region itself which has an assumed cyndrical geometry. All of these c regions are separated by drift regions. c A particle starts its trajectory at the entrance to the first c sextupole (sometimes called 6P1). In four nested loops it is given its C initial velocity and position vectors, while being checked to make sure c its trajectory is allowable by the geometry of the nozzle and the c skimmer. The particle is either focused or defocused due to its c electron spin (M1 in 6P1, M2 in 6P2) within the two sextupoles. c The effects from the nuclear magetic momemt are small, and c disregarding these enables the code to be purely analytic, which c decreases the run time. Thus, this code is ideally c suited for rapid searches to optimize sextupole c placement and geometries. Its cannot determine, for example, the c slight differences for trajectories of hydrogen atoms in states 1 and 2. c After the particle exits the second sextupole it is traced through c the iron plate and into the ionizing region, where hits (particles that c make it to the ionizing region), pathlength (the sum of all pathlengths) c and time in the ionizer as well as a radial distribution (RAD(I)) are c calculated. If a particle exceeds the boundaries of the path geometry c (hits a sextupole, strikes the iron plate, etc.), then the program c loops back to start a new particle through the system. The final c results are weighted with the velocity distributions and output to c the screen. Set P equal to 4 for a file output (see below). c It depends on the characteristics of the ionizer and on one's c application whether total HITS, pathlength(PATH), or time (Q) spent c in the ionizer should be maximized over all initial parameters for best c beam. c The origin of the coordinate system for the analytic equation for c the first sextupole exists where the extrapolated, tapered sextupole c tips would intersect the central beam axis if they were long enough. c c TO RUN THIS PROGRAM: Adjust the parameter list below to your heart's c desire, then just run ATOMSOR. In order to step through a variable, c just insert a DO loop around the program and comment out the variable c in the parameter list. DON'T FORGET to uncomment it when you're done. c If your changes make the run time long, you may want to go to file c output(P=4) and run it as a BATCH job. To find greatest output flux, c maximize HITS. To get best positive ion beam flux out of an ionizer c downstream from your atomic beam source, Pathlenth and Time in the c ionizer(Q) may be the quality factors one wishes to maximize. c c *****Additional notes C WRITEN BY:RON SAYER, JUNE 9,1988 - modified until JUNE 30,1989 C CODE RUNS THROUGH RANGES OF VELOCITY, X-ANGLE, Y-ANGLE, AND RADIUS FOR C A PARTICLE STARTING AT THE ENTRANCE TO THE FIRST SEXTUPOLE C ALL RAYS START ON SEGMENT SO INITIAL X DISPLACEMENT=0 C YI IS INITIAL RADIAL DISTANCE C ALPY IS ANGLE FROM RAY TO Z-AXIS (BEAM AXIS) IN THE Y-Z PLANE C ALPX IS ANGLE FROM RAY TO Z-AXIS (BEAM AXIS) IN THE X-Z PLANE C VZ IS VELOCITY ALONG BEAM AXIS (CONSTANT FOR EACH PARTICLE), C WEIGHTED FROM MEASURED DATA FROM ZURICH, SWITZERLAND C NUC. INST. & METH. A240 (1985) PP.229-236 C ALL OTHER VELOCITIES ARE UNITLESS, SINCE Y AND X ARE FUNCTIONS OF Z, C E.G. VY=DY/DZ=SIN(ALPY) C ----- NOTE: UNITS ARE IN CGS SYSYTEM ------ (YUCK) C C REAL*4 L1,L2,SEP,VZ,THY1,TH2,X,Y,MP,UB REAL*4 ATAPY,LNZS,YI,Y1,VY1,VY2,Y2,BETA REAL*4 FOCUS,IPATHLEN,KINE,INITE,M,BB REAL*4 ALP,IRAD1,ORAD1,RAD2,INC,TNPHASE,PHASE REAL*4 INB1,OUTB1,B2,INCB,INCR,NOZ6PDIS,NOZSKDIS,SKIMRAD,NOZRAD REAL*4 IONPOS,IONRAD,IONLEN,NOZ,NOZX,SKIM,M1,M2 DIMENSION B(2,9),IPATHLEN(9),HIT(9),WEIGHT(9),RAD(20),hita(28) DIMENSION HITR(13) INTEGER STEPMAG,NV,NR,NA,RAYS,P,H C c********************* Parameter List Begins ****************** c P=5 !P=4 FOR FILE OUTPUT, P=5 FOR SCREEN OUTPUT IF (P.EQ.4) THEN OPEN(4,FILE='ATOM.DAT',STATUS='NEW') ENDIF C C ************* CONSTANTS ************* UB=.927E-20 !BOHR MAGNETON (ERG/GAUSS) MP=1.67E-24 !MASS OF PROTON (GRAMS) NUCMASS=1 !1 FOR H AND 2 FOR D C ----SPIN OF PARTICLE M1=+.5 !M=-1/2 MEANS PARTICLES ARE REPELLED FROM Z-AXIS, C !M=+1/2 MEAN PARTICLES ARE FOCUSED in 6P1. M2=+.5 !SAME FOR 6P2 c write(5,*)'enter m1,m2:' c read(5,*)m1,m2 C C ----PRE-SEXTUPOLE REGION DIMENSIONS NOZSKDIS=2.0 !Distance from nozzle to skimmer(cm) NOZ6PDIS=4.0 !Distance from nozzle to 6P1(cm) NOZRAD=0.15 !Nozzle radius(cm) SKIMRAD=0.25 !Skimmmer radius(cm) C C ----SEXTUPOLES PARAMETERS IRAD1=.7 !6P1 entrance pole tip radius from z-axis(cm) ORAD1=1.4 !6P1 exit pole tip radius from z-axis(cm) INB1=7850 !6P1 entrance pole tip field(gauss) write(5,*)'enter INB1:' read(5,*)INB1 C OUTB1=6700 !6P1 exit pole tip field(gauss) B2=6750 !6P2 pole tip field(gauss) RAD2=1.5 !6P2 constant pole tip radius(cm) L1=10.0 !Length of 6P1(cm) L2=10.0 !Length of 6P2(cm) SEP=26.0 !Axial separation of sextupoles(cm) C C ----IONIZER PARAMETERS IONPOS=20.0 !axial distance of iron plate from 6P2(cm) IONLEN=30.0 !effective ionizing region axial length(cm) IONRAD=.75 !radius of ionizing region of the plasma(cm) PLATEHOLE=1.27 !radius of hole in iron plate after 6P2(cm) PLASBEGIN=8.0 !axial distance beyond iron plate where ionizing begins PROBV=85000 !Most probable axial velocity(cm/s) C C ----WEIGHTS FOR WEIGHTED VELOCITY DISTRIBUTION C --MOST PROBABLE VEL. IS AT 4.5 C --EACH INDEX STEPS 10000 CM/S WEIGHT(1)=.05 !35000 CM/S SLOWER THAN PROBV WEIGHT(2)=.267 WEIGHT(3)=.533 WEIGHT(4)=.8 WEIGHT(5)=.8 WEIGHT(6)=.6 WEIGHT(7)=.4 WEIGHT(8)=.2 WEIGHT(9)=.1 !45000 CM/S FASTER THAN PROBV C ----HALF WIDTH DOUBLED -- LEAVE COMMENTED OUT C WEIGHT(1)=.4 !35000 CM/S SLOWER THAN PROBV C WEIGHT(2)=.533 C WEIGHT(3)=.666 C WEIGHT(4)=.8 C WEIGHT(5)=.8 C WEIGHT(6)=.7 C WEIGHT(7)=.6 C WEIGHT(8)=.5 C WEIGHT(9)=.4 !45000 CM/S FASTER THAN PROBV C C C ********************* END OF PARAMETER LIST ****************** C ----ASSIGN Z-AXIS COORDINATES(FOR ANALYTIC FORMULAS) Z1=L1*IRAD1/(ORAD1-IRAD1) !ENTRANCE 6P1 Z2=Z1+L1 !EXIT 6P1 Z3=Z2+SEP !ENTRANCE 6P2 Z4=Z3+L2 !EXIT 6P2 Z5=Z4+IONPOS !IONIZER ENTRANCE Z6=Z5+IONLEN !IONIZER EXIT C C ----INITIALIZE DATA TO ZER0 RAYS=0 !#RAYS SHOT THROUGH ENTRANCE OF 6-POLE 1 C MISSPLAS=0 !RAYS PAST IRON PLATE THAT MISS PLASMA DO 8 I=1,9 !INDEX (I) DENOTES VELOCITY HIT(I)=0 !# PARTICLES THAT GET IN THE IONIZER REGION 8 IPATHLEN(I)=0 !DISTANCE PARTICLE TRAVELS IN IONIZER REGION DO 14 I=1,20 14 RAD(I)=0 !# OF PARTICLES INSIDE RADIUS I*.05CM AT IONIZER ENTR. DO 2050 I=1,28 2050 HITA(I)=0 DO 2060 I=1,13 2060 HITR(I)=0 C C ----FIND CONSTANTS FOR EQNS. OF MOTION (OUTSIDE LOOP--SAVE TIME) DO 7 I=1,9 !FIRST INDEX IS 6-POLE #, SECOND INDEX IS VELOCITY VZ=10000*I+(PROBV-45000) KINE=.5*MP*NUCMASS*VZ*VZ !KINETIC ENERGY OF PARTICLE B(1,I)=UB*INB1*(Z1/IRAD1)*(Z1/IRAD1)/KINE 7 B(2,I)=UB*B2/(RAD2*RAD2*KINE) C LNZS=ALOG(Z2/Z1) !CONSTANT IN EQNS. OF MOTION FOR 6P1 C C ****** MAIN LOOP ****** C ----FOR R (YI),Y-ANGLE (ALPY), AND X-ANGLE (ALPX) C DO 50 YI=0,0.6,.05 !0 TO .60 STEP .05 DO 50 ALPY=-.07,.2,.01 !-.07 TO .2 STEP .01 DO 50 ALPX=0,0.035,.005 !0 TO +.035 STEP .005 C C ----ELIMINATE IMPOSSIBLE ANGLES FOR OUR NOZZLE C AND SKIMMER SET-UP C ----TRACES RAYS BACK TO NOZZLE,USES SMALL ANGLE APPROX. NOZ=ABS(YI-NOZ6PDIS*ALPY) SKIM=ABS(YI-(NOZ6PDIS-NOZSKDIS)*ALPY) NOZX=ABS(NOZ6PDIS*ALPX) SKIMX=ABS((NOZ6PDIS-NOZSKDIS)*ALPX) NOZ=NOZ*NOZ+NOZX*NOZX SKIM=SKIM*SKIM+SKIMX*SKIMX IF (SKIM.GT.(SKIMRAD*SKIMRAD).OR.NOZ.GT.(NOZRAD*NOZRAD)) THEN GO TO 50 ENDIF C C -----INITIALIZE PARAMETERS BEFORE VZ LOOP----- NV=0 !DENOTES VELOCITY NUMBER; I.E. NV=1 C !REFERS TO VZ=PROBV-35000 CM/S (SEE WEIGHTS). VY1=ALPY !SMALL ANGLE APPROXIMATION: EQNS. OF MOTION VX1=ALPX !FOR Y AND X ARE W/RESPECT TO Z; VY1=DY/DZ=TAN(ALPY). C C -----LOOP FOR VELOCITY DISTRIBUTION----------- DO 50 VZ=(PROBV-35000),(PROBV+45000),10000 RAYS=RAYS+1 NV=NV+1 C C ************************ 1ST SEXTUPOLE*************************** C --YI & VY1 ARE Y-POSITION & VELOCITY AT THE ENTRANCE TO 6P1 C --Y2 & VY2 ARE SIMULARLY FOR THE EXIT OF 6P2 C IF (M1.GE.0) THEN C ------------ M=+1/2 -------------- C C -- Y-COMPONENT BETA=SQRT(B(1,NV)-.25) IF (YI.NE.0) THEN TNPHASE=(.5-(Z1*VY1/YI))/BETA PHASE=ATAN(TNPHASE) ATAPY=YI*SQRT((1+TNPHASE*TNPHASE)/Z1) ELSE !TAKE LIMIT AS YI GOES TO 0 PHASE=-1.5708 !PI/2 ATAPY=VY1*SQRT(Z1)/BETA ENDIF THY1=BETA*LNZS+PHASE Y2=ATAPY*SQRT(Z2)*COS(THY1) VY2=Y2/(2*Z2)-ATAPY*BETA*SIN(THY1)/SQRT(Z2) C C -- X-COMPONENT PHASEX=-1.5708 !PI/2 ATAPX=VX1*SQRT(Z1)/BETA THX1=BETA*LNZS+PHASEX X2=ATAPX*SQRT(Z2)*COS(THX1) VX2=X2/(2*Z2)-ATAPX*BETA*SIN(THX1)/SQRT(Z2) C C ----CHECK TO SEE IF RAY HITS 6-POLE TAPER1=(ORAD1-IRAD1)/L1 !TAPER OF 6P1 C ---ANGLE MUST BE GT TAPER TO HIT 6-POLE IF (ALPY.GT.TAPER1. OR. ALPX.GT.TAPER1) THEN C ---CLOSEST APPROACH WHERE SLOPE OF TRAGECTORY=TAPER C ---APPROX. ZMAX NEAR MIDDLE OF 6P1 APPROACH=(TAPER1/ATAPY)*SQRT((Z1+Z2)/(2*B(1,NV))) C ---APPROACH>1 MEANS PARTICLE VERY CLOSE TO 6P1(RARE) C ---GO TO DRIFT IF (APPROACH.GE.1) THEN GO TO 77 !CODE MAY NOW TAKE ACOS(APPROACH) TO FIND SLOPE ENDIF C ACOSAPP=ACOS(APPROACH) PHPRIME=ATAN((B(1,NV)*YI/(Z1*ALPY)-.5)/BETA) ZMAX=Z1*EXP((ACOSAPP-PHPRIME)/BETA) !Z OF CLOSEST APPROACH IF (ZMAX.LT.Z2) THEN !FINDS X AND Y AT ZMAX PHASE2=BETA*ALOG(ZMAX/Z1) THYMAX=PHASE2+PHASE YM=ATAPY*SQRT(ZMAX)*COS(THYMAX) THXMAX=PHASE2+PHASEX XM=ATAPX*SQRT(ZMAX)*COS(THXMAX) ELSE !X AND Y HAVE ALREADY BEEN FOUND AT END XM=X2 ! OF 6P1 YM=Y2 ENDIF R=TAPER1*ZMAX !RADIUS OF 6P1 AT ZMAX IF ((XM*XM+YM*YM).GE.(R*R)) THEN !RAY HITS 6P1 GO TO 50 ENDIF ENDIF C -------------M=-1/2--------------- ELSE BETA=SQRT(B(1,NV)+.25) C --Y COMPONENT BB=(Z1*ALPY+(BETA-.5)*YI)/(2*BETA) CC=BB-YI ZZ=Z2/Z1 Y2=BB*(ZZ)**(BETA+.5)-CC*(ZZ)**(-BETA+.5) C --X COMPONENT BBX=(Z1*ALPX)/(2*BETA) ZZ=Z2/Z1 X2=BBX*(ZZ)**(BETA+.5)-BBX*(ZZ)**(-BETA+.5) C --TEST FOR CRASHING IF ((Y2*Y2+X2*X2).GE.(IRAD1*IRAD1)) THEN GO TO 50 ENDIF VY2=BB*(BETA+.5)*(ZZ)**(BETA-.5) VY2=(VY2+CC*(BETA-.5)*(ZZ)**(-BETA-.5))/Z1 VX2=BBX*(BETA+.5)*(ZZ)**(BETA-.5) VX2=(VX2+CC*(BETA-.5)*(ZZ)**(-BETA-.5))/Z1 ENDIF C ***********************END OF 1ST SEXTUPOLE******************** C C DRIFTING BETWEEN 6-POLES 77 Y3=Y2+VY2*SEP X3=X2+VX2*SEP R=X3*X3+Y3*Y3 !RADIUS OF RAY RR=RAD2*RAD2 IF (R.GE.RR) THEN GO TO 50 !LOST BETWEEN 6-POLES ENDIF C C ************************ 2ND SEXTUPOLE***************************** C --Y3 & VY2(NO VELOCITY CHANGE IN DRIFT SPACE) ARE Y-POS. & VEL. AT C THE ENTRACE TO 6P2: Y4 AND VY4 ARE SIMULARLY Y-POS. & VEL. AT ITS EXIT C IF (M2.GE.0) THEN C ---------------------M=+1/2-------------------- BETA2=SQRT(B(2,NV)) C ----Y-COORDINATE THY2=BETA2*L2 ST=SIN(THY2) CT=COS(THY2) Y4=(VY2/BETA2)*ST+Y3*CT VY4=VY2*CT-BETA2*Y3*ST C ----X-COORDINATE X4=(VX2/BETA2)*ST+X3*CT VX4=VX2*CT-BETA2*X3*ST C C ----CHECK TO SEE IF RAY HITS 6-POLE C --IF PARTICLE STAYS IN "BOUNDS" WITHOUT MAGNET, C --THEN, DON'T HAVE TO CHECK STRAX=X3+VX2*L2 STRAY=Y3+VY2*L2 IF ((STRAY*STRAY+STRAX*STRAX).GT.(RR)) THEN C ---ZMAX IS APPROXIMATION TO CLOSEST APPROACH TO 6-POLE C ---WHERE SLOPE=0 IF (VY2.GE.0) THEN ZMAX=ATAN(VY2/BETA2)/BETA2 !FINDS MAX OR MIN OF ORBIT ELSE ZMAX=(3.141592654+ATAN(VY2/BETA2))/BETA2 ENDIF IF (ZMAX.LT.L2) THEN !FIND X & Y MAX THYMAX2=ZMAX*BETA2 ST=SIN(THYMAX2) CT=COS(THYMAX2) YM=(VY2/BETA2)*ST+Y3*CT XM=VX2*CT-BETA2*X3*ST ELSE !MAX AT END OF 6-POLE XM=X4 YM=Y4 ENDIF IF ((XM*XM+YM*YM).GE.(RR)) THEN GO TO 50 !CRASHES INTO 6P2 ENDIF ENDIF ELSE C -------------------M=-1/2----------------------- C ----Y COORDINATE BBB=SQRT(UB*B2/(RAD2*RAD2*KINE)) A=.5*(Y3+VY2/BBB) C=.5*(Y3-VY2/BBB) Y4=A*EXP(BBB*L2)+C*EXP(-BBB*L2) C ----X COORDINATE BBBX=SQRT(UB*B2/(RAD2*RAD2*KINE)) AX=.5*(X3+VX2/BBBX) CX=.5*(X3-VX2/BBBX) X4=AX*EXP(BBBX*L2)+CX*EXP(-BBBX*L2) C ----CHECK FOR CRASHING IF ((X4*X4+Y4*Y4).GE.(RAD2*RAD2)) THEN GO TO 50 ENDIF VY4=A*BBB*EXP(BBB*L2)-C*BBB*EXP(-BBB*L2) VX4=AX*BBBX*EXP(BBBX*L2)-CX*BBBX*EXP(-BBBX*L2) ENDIF C C -----JUST IN CASE V(PURP)=0 THEN DRIFT IF (VY4*VY4+VX4*VX4.EQ.0) THEN IF (X4*X4+Y4*Y4.GT.IONRAD*IONRAD) THEN GO TO 50 ELSE HIT(NV)=HIT(NV)+1 I=ALPY*100+8.01 HITA(I)=HITA(I)+1 I=YI*20+1.01 HITR(I)=HITR(I)+1 IPATHLEN(NV)=IPATHLEN(NV)+IONLEN GO TO 50 ENDIF ENDIF C ************************ END OF 2ND SEXTUPOLE ******************** C C ----DRIFTS TO ENTRANCE OF IONIZER Y5=Y4+VY4*IONPOS Y6=Y5+VY4*PLASBEGIN X5=X4+VX4*IONPOS X6=X5+VX4*PLASBEGIN R5=SQRT(X5*X5+Y5*Y5) !RADIUS AT IRON PLATE AFTER 6P2 R6=SQRT(X6*X6+Y6*Y6) !RADIUS AT BEGINING OF PLASMA C C ***************************** IONIZER ************************ C ----FIND RADIAL DENSITY---- IR=(20*R5)+1 IF (IR.LE.20) THEN RAD(IR)=RAD(IR)+WEIGHT(NV) ENDIF C C ----REJECTS IF PARTICLE HITS IRON PLATE OR MISS PLASMA---- IF ((ABS(R5).LT.PLATEHOLE).AND.(ABS(R6).GT.IONRAD)) THEN MISSPLAS=MISSPLAS+1 GO TO 50 ENDIF IF ((ABS(R5).GE.PLATEHOLE).OR.(ABS(R6).GT.IONRAD)) THEN GO TO 50 ENDIF C C ----FINDS PATHLENGTH IN IONIZER---- C --FINDS PARTICLES INTERSECTION WITH IONIZING CYLINDER RR=IONRAD*IONRAD IF (VX4.NE.0) THEN M=VY4/VX4 !SLOPE BB=Y6-M*X6 !Y-INTECEPT IF (VX4.GE.0) THEN X=(-BB*M+SQRT(RR*(M*M+1)-BB*BB))/(M*M+1) ELSE X=(-BB*M-SQRT(RR*(M*M+1)-BB*BB))/(M*M+1) ENDIF Z=(X-X6)/VX4 ELSE X=X6 Y=(ABS(VY4)/VY4)*SQRT(RR-X*X) Z=(Y-Y6)/VY4 ENDIF C C ----ADDS UP PATHLENGTH:PLENGTH (APPROXIMATELY)=SQRT(Z*Z+X*X+Y*Y) IF (Z.GE.IONLEN) THEN IPATHLEN(NV)=IPATHLEN(NV)+IONLEN ELSE Y=Y6+VY4*PLENGTH XD=X-X6 YD=Y-Y6 IPATHLEN(NV)=IPATHLEN(NV)+PLENGTH ENDIF HIT(NV)=HIT(NV)+1 !COUNTS HITS INTO IONIZER I=ALPY*100+8.01 HITA(I)=HITA(I)+1 I=YI*20+1.01 HITR(I)=HITR(I)+1 C *********************** END OF IONIZER ********************** C 50 CONTINUE !END OF LOOP FOR VZ,ALPX,ALPY AND YI C C ---INITIALIZE OUTPUT TO ZERO Q=0 PATH=0 H=0 WH=0 WPATH=0 WQ=0 C C ---ADD UP UNWEIGHTED PARAMETERS AND ADD WEIGHT DO 244 I=1,9 VZ=PROBV+(I-4.5)*10000 H=H+HIT(I) PATH=PATH+IPATHLEN(I) Q=Q+IPATHLEN(I)/VZ WH=WH+HIT(I)*WEIGHT(I) WPATH=WPATH+IPATHLEN(I)*WEIGHT(I) 244 WQ=WQ+IPATHLEN(I)*WEIGHT(I)/VZ C C ---PRINT RESULTS TO SCREEN (P=5) OR FILE (P=4) C Q=Q/RAYS !OPTIONAL FOR QUALITY VALUE C c********** OUTPUT ********************** c WRITE(P,*)'DOING HITS IN Vel. 500-1300m/s' DO 2000 I=1,9 2000 WRITE(P,*)HIT(I) C WRITE(P,*)'# RAYS= ',RAYS WRITE(P,*)'IONIZER POSITION=',(IONPOS+PLASBEGIN) C WRITE(P,*)'B1 IN = ',INB1 C WRITE(P,*)'B2 = ',B2 WRITE(P,*)'MOST PROBABLE VELOCITY= ',(PROBV/100),' M/S' C WRITE(P,*)'LENGTH OF 6-POLE #1= ',L1 C WRITE(P,*)'# HITS= ',H C WRITE(P,*)'TOTAL PATHLENGTH = ',PATH C WRITE(P,*)'QUALITY FACTOR = ',Q C WRITE(P,*)'Q IS TIME IN IONIZER/#RAYS SHOT THROUGH' C WRITE(P,*)'# MISSING PLASMA = ',MISSPLAS WRITE(P,*)'WEIGHTED PARAMETERS:' WRITE(P,*)' HITS: ',WH C WRITE(P,*)' PATHLENGTH: ',WPATH C WRITE(P,*)' Q: ',WQ C WRITE(P,*)'WEIGHTED HITS ',IONPOS,' AFTER 6P1:' 4444 WRITE(P,*)' ' c IF (1.LT.2) THEN c GO TO 4444 c ENDIF IF (P.EQ.4) THEN CLOSE(4) ENDIF STOP END