C
C
C				A T O M S T A T E
C
C
C                ********** INTRODUCTION ****************
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 moment are usually small, unless
c	the particle is in a low magnetic field (near the central beam axis).
c	The magnetic moment of the atom can then change for certain spin states
c	of hydrogen and deuterium.  These changes in states 2 and 4 for
c	Hydrogen(2,3,5,6 FOR D) are taken into account by stepping through
c	6P1 and 6P2 STEPMAG # of times, but there is significantly
c	longer runtime than if the code ignored this effect.
c	If this is a small effect for the user, or one only is running states
c	1 & 3 for H(1 & 4 for D) then STEPMAG can be set to 1 and the program
c	will run very quickly, since an analytic formula is programmed in for
c	both 6P1 (a tapered sextupole) and 6P2.  If this is done, CRASHING
c	CONDITIONS MUST be changed in 6P1 and 6P2 (particle could crash into
c	pole-tips and still come out inside the output radius of a sextupole
c	in this code if the atom's position isn't checked often enough(STEPMAG
c	too small)).
c	  NOTE:  If it is the user's intent to optimize one's atomic beam
c	system and the user suspects that the changing magnetic moment of
c	hydrogen and dueterium at low magnetic fields isn't an overriding
c	consern, I SUGGEST that the user run ATOMSOR, a code that doesn't
c	take into account the changing magnetic moment, yet RUNS FASTER
c	than this code.  To calculate polarizations, however, it is nessessary
c	to use this code.
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 TUMAG.  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	U(MAGNETIC MOMENT) IS NO LONGER CONSTANT!!!!!!
c
C	WRITEN BY:RON SAYER-FUN AND ZZZZ MAY 19,1989 - last modifed Sept.18,1989
C	CODE RUNS THROUGH RANGES OF VELOCITY, ALPHA, AND RADIUS FOR A 
C	  PARTICLE STARTING AT THE ENTRANCE TO THE 1rst 6-POLE.
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)
C	L1 & L2 ARE LENGTHS OF 1rst and 2nd MAGNETS (SEP IS MAGNET SEPARATION)
C	LENGTH IS DISTACE FROM ENTRANCE OF 6-POLE1 TO EXIT OF 6-POLE2
C	MISSPLAS IS # OF PARTICLES THAT GO THROUGH IRON PLATE BUT MISS 
C		THE IONIZING REGION AT PLASBEGIN
C	RAD STORES HITS INSIDE RADI OF .05,.1,...,1.0
C	VZ=VELOCITY(ALONG BEAM AXIS)=550 TO 1350 STEP 100
C	USED MEASURED DATA FROM ZURICH,SWITZERLAND FOR WEIGHTS
C		NUC. INST. & METH.  A240 (1985) pp.229-236
C	ALL OTHER VELOCITIES ARE UNITLESS (E.G. VY1=VEL. IN Y/ VZ)
C	----NOTE ALL UNITS ARE IN CGS SYSTEM (YUCK)
C
	REAL*4 L1,L2,SEP,VZ,TH1,TH2,X,Y,MP,UB
	REAL*4 FOCUS,IPATHLEN,INITE,M,BB
	REAL*4 ALP,IRAD1,ORAD1,RAD2,INC,TNPHASE,PHASE
	REAL*4 INB1,OUTB1,B2,NOZ6PDIS,NOZSKDIS
	REAL*4 IONPOS,IONRAD,IONLEN,NOZ,NOZX,SKIM,M1,M2,NOZRAD,SKIMRAD
	REAL*4 NP,N0,NM
	DIMENSION B(2,9),IPATHLEN(9),HIT(9),WEIGHT(9),RAD(20),RKINE(9)
	DIMENSION WHSTATE(6,2),CURRENT(7),PVECTOR(7),PTENSOR(7)
	DIMENSION PV(14,7),PT(14,7),CUR(14,7)
	INTEGER STEPMAG,NV,NR,NA,RAYS,P,H,STATE,BT1TO4,BT3TO5
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='DSTATE.DAT',STATUS='NEW')
	  ENDIF
C
C		***********CONSTANTS*********
	UB=.927E-20	!BOHR MAGNETON (ERG/GAUSS)
	MP=1.67E-24	!MASS OF PROTON (GRAMS)
C
C		----SPIN OF PARTICLE
C		----PRE SEXTUPOLE 1 DIMENSIONS
	NOZ6PDIS=4.	!Distance from nozzle to skimmer (cm) 
	NOZSKDIS=2.	!Distance from 6P1 (cm)
	SKIMRAD=0.25	!Skimmer radius (cm)
	NOZRAD=0.15	!Nozzle radius (cm)
C
C		----IONIZER PARAMETERS
	IONPOS=20.0	!axial distance of iron plate from 6P2 exit(cm)
	IONLEN=30.	!effective ionizing region length(cm)
	IONRAD=.75	!radius of IONIZING PART OF THE PLASMA(cm)
	PLATEHOLE=1.27	!RADIUS OF HOLE IN IRON PLATE AFTER 6P2 (cm)
	PLASBEGIN=8.0	!Axial distance beyond iron WHERE PLASMA Begins(cm)
C
	PROBV=60000	!CM/S
C
c		----SEXTUPOLE PARAMETERS
	IRAD1=.7	!6P1 entrance pole tip radius in CM
	ORAD1=1.4	!6P1 exit pole tip radius in CM
	INB1=7850.	!6P1 entrance pole tip field in gauss
	OUTB1=6700.	!6P1 exit pole tip field in gauss
	B2=6750.	!6P2 pole tip field in Gauss
	RAD2=1.5	!6P2 constant pole tip radius
	L1=10.		!LENGTH OF 6P1 (cm)
	L2=10.		!length of 6P2 (cm)
	SEP=26.		!axial SEPARATION OF SEXTUPOLES(cm)
	STEPMAG=10	!#OF STEPS THROUGH SEXTUPOLE(FOR CHANGING U)
C
C		--NUCLEAR MASS AND STATE
	NUCMASS=2	!NUCLEAR MASS(1 FOR H AND 2 FOR D)
	IF (NUCMASS.EQ.1) THEN
	  Bo=507.	!hydrogen's b not(for changing magnetic moment)
	ELSE
	  Bo=117.	!deuterium's b not(for changing magnetic moment)
	  ENDIF
	STATE=1
c		--DEUTERIUM'S TRANSITIONS (0 IS ON AND 1 IS OFF)
	BT1TO4=0	!TRANSITION STATE 1 TO 4 BETWEEN SEXTUPOLES
	BT3TO5=0	!TRANSITION STATE 3 TO 5 BETWEEN SEXTUPOLES
c	--Statements assign values to M1 and M2 from State
c	    M1 is the Electron spin's direction(+ for up, - for down) in 6P1
c	    M2 is the Electron spin's direction(+ for up, - for down) in 6P2
c
	IF (NUCMASS.EQ.1) THEN
C	---- HYDROGEN ----
	  IF (STATE.LE.2) THEN
	    M1=1
	    M2=1
	  ELSE
	    M1=-1
	    M2=-1
	    ENDIF
	ELSE
C	---- DEUTERIUM ----
	  IF (STATE.LE.3) THEN
	    M1=1
	    IF ((STATE*BT1TO4).EQ.1.OR.(STATE*BT3TO5).EQ.3) THEN
	      M2=-1
	    ELSE
	      M2=1
	      ENDIF
	  ELSE
	    M1=-1
	    IF ((STATE*BT1TO4).EQ.4.OR.(STATE*BT3TO5).EQ.5) THEN
	      M2=1
	    ELSE
	      M2=-1
	      ENDIF
	    ENDIF
	  ENDIF	
C
C	--WEIGHTS FOR WEIGHTED VELOCITY DISTRIBUTION
C		--MOST PROBABLE VEL. IS AT 4.5
C		--EACH INDEX STEPS 100 M/S
	WEIGHT(1)=.05		!350 M/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		!450 M/S FASTER THAN PROBV
C	----HALF WIDTH DOUBLED
C	WEIGHT(1)=.4		!350 M/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		!450 M/S FASTER THAN PROBV
c
c***************** End Of Parameter List ****************
c
C	   ----ASSIGN Z-AXIS COORDINATES (FOR ANALYTIC FORMUALAS)
	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	--THIS LOOP CALCUALTES KINETIC ENERGY & B (SAVES RUN TIME)
	DO 888 I=1,9		!RKINE IS KINETIC ENERGY
	  VZ=PROBV+(I-4.5)*10000
	  RKINE(I)=.5*MP*NUCMASS*VZ*VZ
888	  B(2,I)=UB*B2/(RAD2*RAD2*RKINE(I))

C
	RINC=(ORAD1-IRAD1)/STEPMAG	!STEPS FOR RADIUS IN 6P1
	BINC=(OUTB1-INB1)/STEPMAG	!STEPS FOR B-FIELD IN 6P1
	ZINC1=L1/STEPMAG		!STEPS FOR AXIAL POSITION IN 6P1
	ZINC2=L2/STEPMAG		!STEPS FOR AXIAL POSITION IN 6P2
C
C    *********INITIALIZE DATA TO ZERO*******
	RAYS=0		!#RAYS SHOT THROUGH ENTRANCE OF 6-POLE 1
	MISSPLAS=0
C
C	----EACH INDEX COORESPONDS TO A VELOCITY WHICH LATER IS WEIGHTED
	DO 8 I=1,9
	  HIT(I)=0		!# PARTICLES THAT ENTER IONIZING REGION
8	  IPATHLEN(I)=0		!DISTANCE PARTICLE TRAVELS IN IONIZER
C
	DO 14 I=1,20
14	  RAD(I)=0		!# PARTICLES INSIDE I*.05 AT IONIZER ENTRANCE
C
C	MAIN LOOP !!!!!!!!!!!!!!!!!!! FOR YI,ALPY,ALPX, AND VZ
C
	DO 50 YI=0.0,.6,.05		!0 TO .60 STEP .05
C	WRITE(5,*)YI
	  DO 50 ALPY=-.07,.2,.01	!-.07 TO .2 STEP .01
	  DO 50 ALPX=0.0,.035,.005	!0    TO +.035 STEP .005
C		----ELIMINATES IMPOSSIBLE ANGLES FOR OUR NOZZLE
C	        ----     AND SKIMMER SET-UP
C		----TRACES RAYS BACK TO NOZZLE,USES SMALL ANGLE APPROX.
C		----DISTANCES:NOZZLE TO SKIMMER = 2.0 CM
C		----DISTANCES:NOZZLE TO 6-POLE1=NOZ6PDIS
	    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
C		--IF TRACED RAY CANNOT HAVE COME FROM NOZZLE DUE TO
C		--     NOZZLE & SKIMMER GEOMETRY, THEN GOTO NEXT RAY
	    IF (SKIM.GT.(SKIMRAD*SKIMRAD).OR.NOZ.GT.(NOZRAD*NOZRAD)) THEN
	      GO TO 50
	      ENDIF
C
	    NV=0
C
C	------LOOP FOR VELOCITY DISTRIBUTION-----
	    DO 50 VZ=(PROBV-35000),(PROBV+45000),10000
	      RAYS=RAYS+1
	      NV=NV+1
C
C	******************** 1rst SEXTUPOLE ********************
C	--SET INITIAL CONDITIONS FOR LOOP
	      VY1=ALPY		!VY1 ACTUALLY = SIN(ALPY) BUT SMALL
	      VX1=ALPX		!APPROXIMATION GOOD ENOUGH
	      XI=0
	      X2=0
	      Y2=YI		
	      R=IRAD1		!RADIUS OF SEXTUPOLE
	      B1=INB1		!B-FIELD OF SEXTUPOLE AT PRESENT AXIAL POSITION
	      YII=YI	!Y-POSTION		--(FOR LINEAR B-FIELD)
	      PRADSQ=X2*X2+Y2*Y2	!RADIUS OF PARTICLE
	      RSQ=R*R			!RADIUS OF SEXTUPOLE
	      IF (M1.GE.0) THEN
C	--------------- M=+1/2 ----------------
C	** AXIALLY STEPS DOWN THE SEXTUPOLE **
C
		DO 400 Z=Z1,(Z2-ZINC1*.5),ZINC1
C
		SQZZ=SQRT(Z+ZINC1)	!DUMMY VARIABLE FOR SAVING RUNNING TIME
		ALLL=ALOG(1+(ZINC1/Z))
	  	B(1,NV)=UB*B1*(Z/R)*(Z/R)/RKINE(NV)
C
C	-----ACCOUNTS FOR NON-CONSTANT u
C	--MAGNETIC MOMENT OF ATOM DEPENDS ON PARTICLES RADIAL POSITION
		IF (NUCMASS.EQ.1) THEN
		  IF (STATE.EQ.2) THEN
		    FRAC=(B1/Bo)*(PRADSQ/RSQ)
		    B(1,NV)=B(1,NV)*FRAC/SQRT(FRAC*FRAC+1)
		    ENDIF
		ELSE
		  IF (STATE.EQ.2) THEN
		    FRAC=(B1/Bo)*(PRADSQ/RSQ)
		    B(1,NV)=B(1,NV)*(FRAC+.33333)/SQRT(FRAC*FRAC+FRAC*.66666+1)
		    ENDIF
		  IF (STATE.EQ.3) THEN
		    FRAC=(B1/Bo)*(PRADSQ/RSQ)
		    B(1,NV)=B(1,NV)*(FRAC-.33333)/SQRT(FRAC*FRAC-FRAC*.66666+1)
		    ENDIF
		  ENDIF
C
C	--IN CASE OF VERY WEAK FIELD, LET PARTICLE DRIFT
		IF ((B(1,NV)-.251).GE.0) THEN 
 	          BETA=SQRT(B(1,NV)-.25)
		ELSE
		  Y2=YII+VY2*ZINC1
		  X2=XI+VX2*ZINC1
		  GOTO 900
		  ENDIF
C
C	---Y-COMPONENT	
	        IF (YII.NE.0) THEN
	          TNPHASE=(.5-(Z*VY1/YII))/BETA
	          PHASE=ATAN(TNPHASE)
	          A=YII*SQRT((1+TNPHASE*TNPHASE)/Z)
	        ELSE			!LET YI GOTO 0
		  PHASE=-1.5708	!PI/2
		  A=VY1*SQRT(Z)/BETA
		  ENDIF
	        TH1=BETA*ALLL+PHASE
	        Y2=A*SQZZ*COS(TH1)
	        VY2=Y2/(2*(Z+ZINC1))-A*BETA*SIN(TH1)/SQZZ
C	  X-COMPONENT
		IF (XI.NE.0) THEN
		  TNPHASEX=(.5-(Z*VY1/XI))/BETA
		  PHASEX=ATAN(TNPHASEX)
		  AX=XI*SQRT((1+TNPHASEX*TNPHASEX)/Z)
		ELSE
	          PHASEX=-1.5708	!PI/2
	          AX=VX1*SQRT(Z)/BETA
		  ENDIF
	        THX1=BETA*ALLL+PHASEX
	        X2=AX*SQZZ*COS(THX1)
	        VX2=X2/(2*(Z+ZINC1))-AX*BETA*SIN(THX1)/SQZZ
C
C		----CHECK TO SEE IF RAY HITS 6-POLE
900		R=R+RINC
		B1=B1+BINC
		PRADSQ=X2*X2+Y2*Y2	!RADIUS OF PARTICLE
		RSQ=R*R			!RADIUS OF SEXTUPOLE
		IF (PRADSQ.GE.RSQ) THEN
	  	  GO TO 50
		  ENDIF
		YII=Y2
		XI=X2
		VX1=VX2
		VY1=VY2
400 		CONTINUE
C	---------------- M=-1/2 ----------------------
		ELSE
		  DO 410 Z=Z1,(Z2-ZINC1*.5),ZINC1
C
	  	  B(1,NV)=UB*B1*(Z/R)*(Z/R)/RKINE(NV)
C
C	-----ACCOUNTS FOR NON-CONSTANT u
C	--MAGNETIC MOMENT OF ATOM DEPENDS ON PARTICLES RADIAL POSITION
		IF (NUCMASS.EQ.1) THEN
		  IF (STATE.EQ.4) THEN
		    FRAC=(B1/Bo)*(PRADSQ/RSQ)
		    B(1,NV)=B(1,NV)*FRAC/SQRT(FRAC*FRAC+1)
		    ENDIF
		ELSE
		  IF (STATE.EQ.5) THEN
		    FRAC=(B1/Bo)*(PRADSQ/RSQ)
		    B(1,NV)=B(1,NV)*(FRAC-.33333)/SQRT(FRAC*FRAC-FRAC*.66666+1)
		    ENDIF
		  IF (STATE.EQ.6) THEN
		    FRAC=(B1/Bo)*(PRADSQ/RSQ)
		    B(1,NV)=B(1,NV)*(FRAC+.33333)/SQRT(FRAC*FRAC+FRAC*.66666+1)
		    ENDIF
		  ENDIF
C
		  IF ((B(1,NV)+.25).LE.0) THEN
		    Y2=YII+VY2*ZINC1
		    X2=XI+VX2*ZINC1
		    GOTO 905
		  ELSE
        	    BETA=SQRT(B(1,NV)+.25)
		    ENDIF
C
C	***** THESE 5 STEPS JUST SAVE RUN TIME FOR THE FORMULAS BELOW ******
		  BETAMP5=BETA-.5
		  BETAPP5=BETA+.5
        	  ZZ=1+ZINC1/Z
		  BPTP=(ZZ)**BETAPP5
		  BMTP=(ZZ)**BETAMP5
C
C	  -Y COMPONENT
         	  BB=(Z*VY1+BETAMP5*YII)/(2*BETA)
        	  CC=BB-YII
         	  Y2=BB*BPTP-CC/BMTP
        	  VY2=(BB*BETAPP5*BMTP+CC*BETAMP5/BPTP)/Z
C	  -X COMPONENT
         	  BBX=(Z*VX1+BETAMP5*XI)/(2*BETA)
		  CCX=BBX-XI
         	  X2=BBX*BPTP-CCX/BMTP
        	  VX2=(BBX*BETAPP5*BMTP+CCX*BETAMP5/BPTP)/Z
C	--TEST FOR CRASHING
905		  R=R+RINC
		  B1=B1+BINC
		  PRADSQ=X2*X2+Y2*Y2
		  RSQ=R*R
          	  IF (PRADSQ.GE.RSQ) THEN
            	    GO TO 50
            	    ENDIF
		  YII=Y2
		  XI=X2
		  VX1=VX2
410		  VY1=VY2
		  ENDIF
C**************** END OF FIRST SEXTUPOLE *****************
C
C	DRIFTING BETWEEN 6-POLES
77	      Y3=Y2+VY2*SEP
	      X3=X2+VX2*SEP
	      R=X3*X3+Y3*Y3	!RADIUS^2 OF RAY
	      RR=RAD2*RAD2
	      IF (R.GE.RR) THEN
	        GO TO 50	!LOST BETWEEN 6-POLES
	        ENDIF
	      X4=X3
	      Y4=Y3
	      PRADSQ=X4*X4+Y4*Y4
C
C   ********************* 2nd 6-POLE **********************
C	--Y3 AND VY2(NO VELOCITY CHANGE IN DRIFT SPACE) ARE Y-POS. & VEL. AT
C	THE ENTRANCE TO 6P2:  Y4 AND VY4 ARE SIMULARLY Y-POS. & VEL. AT ITS EXIT
C
	      IF (M2.GE.0) THEN
C	---------------- M=+1/2 -------------------
		DO 420 Z=0,(L2-ZINC2*.5),ZINC2
C
		BBB=B(2,NV)
C	-----ACCOUNTS FOR NON-CONSTANT u
		IF (NUCMASS.EQ.1) THEN
		  IF (STATE.EQ.2) THEN
		    FRAC=B2/Bo*(PRADSQ/(RAD2*RAD2))
		    BBB=BBB*FRAC/SQRT(FRAC*FRAC+1)
		    ENDIF
		ELSE
		  IF (STATE.EQ.2) THEN
		    FRAC=B2/Bo*(PRADSQ/(RAD2*RAD2))
		    BBB=BBB*(FRAC+.33333)/SQRT(FRAC*FRAC+FRAC*.66666+1)
		    ENDIF
		  IF (STATE.EQ.3) THEN
		    FRAC=B2/Bo*(PRADSQ/(RAD2*RAD2))
		    BBB=BBB*(FRAC-.33333)/SQRT(FRAC*FRAC-FRAC*.66666+1)
		    ENDIF
		  ENDIF
C
		IF (BBB.LE.0) THEN	!DRIFT
		  Y4=Y3+VY2*ZINC2
		  X4=X3+VX2*ZINC2
		  GOTO 910
		  ENDIF
	        SB=SQRT(BBB)
C			Y-COORDINATE
 	        TH2=SB*ZINC2
	        ST=SIN(TH2)
	        CT=COS(TH2)
	        Y4=(VY2/SB)*ST+Y3*CT
	        VY4=VY2*CT-SB*Y3*ST
C			X-COORDINATE
	        X4=(VX2/SB)*ST+X3*CT
	        VX4=VX2*CT-SB*X3*ST
C
910		PRADSQ=X4*X4+Y4*Y4
C		----CHECK TO SEE IF RAY HITS 6-POLE
	          IF (PRADSQ.GE.(RR)) THEN
	    	    GO TO 50	!CRASHES INTO 6P2
		    ENDIF
		Y3=Y4
		X3=X4
		VX2=VX4
420		VY2=VY4
	      ELSE
C	-------------- M=-1/2 ------------
		DO 430 Z=0,(L2-ZINC2*.5),ZINC2
C
		BBB=B(2,NV)
C
C	-----ACCOUNTS FOR NON-CONSTANT u
		IF (NUCMASS.EQ.1) THEN
		  IF (STATE.EQ.4) THEN
		    FRAC=B2/Bo*(PRADSQ/(RAD2*RAD2))
		    BBB=BBB*FRAC/SQRT(FRAC*FRAC+1)
		    ENDIF
		ELSE
		  IF (STATE.EQ.5) THEN
		    FRAC=B2/Bo*(PRADSQ/(RAD2*RAD2))
		    BBB=BBB*(FRAC-.33333)/SQRT(FRAC*FRAC-FRAC*.66666+1)
		    ENDIF
		  IF (STATE.EQ.6) THEN
		    FRAC=B2/Bo*(PRADSQ/(RAD2*RAD2))
		    BBB=BBB*(FRAC+.33333)/SQRT(FRAC*FRAC+FRAC*.66666+1)
		    ENDIF
		  ENDIF
C
		IF (BBB.LE.0) THEN
		  Y4=Y3+VY2*ZINC2
		  X4=X3+VX2*ZINC2
		  GOTO 920
		  ENDIF
      		BBB=SQRT(BBB)
		EBZ=EXP(BBB*ZINC2)
C			--Y COORDINATE
       		A=.5*(Y3+VY2/BBB)
       		C=.5*(Y3-VY2/BBB)
       		Y4=A*EBZ+C/EBZ
     		VY4=A*BBB*EBZ-C*BBB/EBZ
C			--X COORDINATE
        	AX=.5*(X3+VX2/BBB)
       		CX=.5*(X3-VX2/BBB)
       		X4=AX*EBZ+CX/EBZ
     		VX4=AX*BBB*EBZ-CX*BBB/EBZ
C			----CHECK FOR CRASHING
920		PRADSQ=X4*X4+Y4*Y4
          	IF (PRADSQ.GE.RR) THEN
            	  GO TO 50
            	  ENDIF
		Y3=Y4
		X3=X4
		VX2=VX4
430		VY2=VY4
		ENDIF
C****************** END OF SECOND SEXTUPOLE ***********************
C
C		-----JUST IN CASE V(PURP)=0 THEN JUST 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
		  IPATHLEN(NV)=IPATHLEN(NV)+IONLEN
		  GO TO 50
		  ENDIF
		ENDIF
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	----FOR RADIAL DENSITY----
C	** PUT PARTICLES IN THEIR RESPECTIVE ANULUSES **
	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--
C	** ATOMS THAT GO THROUGH THE IRON PLATE BUT MISS THE PLASMA **
	      IF ((R5.LT.PLATEHOLE).AND.(R6.GT.IONRAD)) THEN
		MISSPLAS=MISSPLAS+1
	        ENDIF
C	** ATOMS THAT HIT THE IRON PLATE OR MISS THE PLASMA **
	      IF ((R5.GE.PLATEHOLE).OR.(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:Z (APPROXIMATELY)=SQRT(Z*Z+X*X+Y*Y)
	      IF (Z.GE.IONLEN) THEN
		IPATHLEN(NV)=IPATHLEN(NV)+IONLEN
	      ELSE
		Y=Y6+VY4*Z
		XD=X-X6
		YD=Y-Y6
		IPATHLEN(NV)=IPATHLEN(NV)+Z
		ENDIF
	      HIT(NV)=HIT(NV)+1		!COUNTS HITS INTO IONIZER
50	  CONTINUE
C	********************* END OF IONIZER *******************
C
C	--INITIALIZE OUTPUT VARIABLES TO ZERO
	Q=0
	PATH=0
	H=0
	WH=0
	WPATH=0
	WQ=0
C		--ADDING UP UNWEIGHTED HITS AND PATHLENGTH
	DO 244 I=1,9
	  VZ=PROBV+(I-4.5)*10000	!VZ IS VELOCITY IN Z-DIRECTION
	  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	Q=Q/RAYS
C
C************************* OUTPUT ************************
C
	WRITE(P,*)'STATE= ',STATE
	WRITE(P,*)'TRANS 1 T 4  &  3 T 5= ',BT1TO4,BT3TO5
	WRITE(P,*)'STEPMAG= ',STEPMAG
	WRITE(P,*)'# RAYS= ',RAYS
C	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
	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
	WRITE(P,*)'  PATHLENGTH: ',WPATH
	WRITE(P,*)'  Q: ',WQ
C	WRITE(P,*)'WEIGHTED HITS ',IONPOS,' AFTER 6P1:'
4444	WHSTATE(STATE,NTRANS)=WQ*100
	WRITE(P,*)' '
c
	IF (P.EQ.4) THEN
	  CLOSE(4)
	  ENDIF
	STOP
	END