c------------------------------------------------------------------------
c	Raytracing of 3-d motion in polarized target field by solution
c	of differential equations of motion via 4th order Runge-Kutta. Field 
c	orientation arbitrary. Horizontal plane is (x,y).
c
c	Note: time is measured in nanoseconds through-out.
c------------------------------------------------------------------------
c
!	beam along x -OARA 12/97
!	modified to allow electrons along beam
!	and to select particles only from origin, not random
!	Also, double precision variables, some constants
!	(such as arguments passed to functions - no /autodbl)

	implicit real*8(a-h,o-z)
	dimension y(9),yout(9),dydx(9),v_i(3),w_f(3),v_i_save(3)
	dimension B_field_z(51,51),B_field_r(51,51),zz(51),rr(51)
	dimension prcs(3),prco(3)
	dimension B_axis(3),B_lab(3)
	dimension BETA_axis(3),BETA_lab(3)
	logical*1 uniform_test,recoils,backwards,back_comp
	logical*1 view_x,view_y,view_z
	character*1 ans,forb,bnom,ansp,ansq,anse,q_e,q_sign
	real*8 K
	dimension ywatch(9),dywatch(9)
c
	common/field_strength/B_field_z,B_field_r,zz,rr
	common/field_angles/tht_b,phi_b
	common/field_type/uniform_test
	common/conversion_factor/factor
	common/view/view_x,view_y,view_z
        common/back/backwards,back_comp
c
	data tgt_r/1.25/			!radius of tgt. material
	data iseed/-12311/			!for random number generator
	data g22/0.001159652d0/
	external sind,cosd,tand,asind,acosd,atand
c
	pi       = 3.141592654d0
	c        = 30.d0		!speed of light (cm/nanosec)
	pmass    = 938.28d0		!proton mass in MeV
	dmass    = 1875.6d0           !deuteron mass in MeV
	emass    = 0.511d0		!electron   "      "
c
	uniform_test = .false.
	recoils      = .true.
	backwards    = .false.
	view_x       = .false.
	view_y       = .false.
	view_z       = .false.
	B_max        = 5.1d0
	back_comp    = .false.
c
c	Initial spin direction
c
	do i = 1,3
	prcs(i) = 0.d0
	enddo
c
c	Ask questions
c
	write(6,'(a)')' Enter E_0 (GeV), theta_e and N_events >'
	read(5,*)e0_e,tht_e,n_events
	write(6,'(a)')
	1	' Enter t_max (ns), # of time steps and phi_step >'
	read(5,*)t_max,n_steps,phi_step
	write(6,'(a)')' Enter spherical angles of B (tht,phi) >'
	read(5,*)tht_b,phi_b
	write(6,'(a)')' Random start point (y) or from origin (n) >'
	read(5,'(a)')ans
	origin = 0.d0
	if(ans .eq. 'y' .or. ans .eq. 'Y') origin = 1.0d0
	write(6,'(a)')' Uniform field test? (y/n) >'
	read(5,'(a)')ans
	if(ans .eq. 'y' .or. ans .eq. 'Y')uniform_test = .true.
	if(.not. uniform_test) then
	 write(6,'(a)')' Nominal central field (5.1 T) (y/n) >'
	 read(5,'(a)')bnom
	 if(bnom .ne. 'y' .and. bnom .ne. 'Y') then
	  write(6,'(a)')' Central field [T] >'
	  read(5,*) B_max
	 end if
	end if
	write(6,'(a)')' Track electrons or recoil particles? (e/r) >'
	read(5,'(a)')anse
	if(anse .eq. 'e' .or. anse .eq. 'E') then
	 recoils = .false.
	 write(6,'(a)')' Electron charge? (n[eg]/p[os]) >'
	 read(5,'(a)')q_e
	 charge = 1.d0
	 if(q_e .eq. 'p' .or. q_e .eq. 'P') charge = -1.d0
	endif
	write(6,'(a)')
	1	' Is target proton, deuteron or electron? (p/d/e) >'
	read(5,'(a)')ansp
	if(ansp .eq. 'p' .or. ansp .eq. 'P')then
	 tmass = pmass
	 write(6,'(a)')
	1	' Proton q-elastic (in 2H) or inelastic (free p)? (q/i) >'
	 read(5,'(a)')ansq
	 if(ansq .eq. 'q' .or. ansp .eq. 'Q')then
	  ebind  = 2.223				! e-p q. el.
	 else
	  write(6,'(a)')' Enter invariant mass W [GeV] (Mp OK too) >'
	  read(5,*)W
	  ebind  = 0.d0					! e-d elast.
	 end if
	else if(ansp .eq. 'd' .or. ansp .eq. 'D')then
	 tmass = dmass
	 ebind  = 0.d0					! e-d elast.
	else
	 tmass = emass					! Moeller
	 ebind  = 0.d0
	 ansp = 'e'
	endif
	if(.not. recoils)then
	     write(6,'(a)')' Enter radial offset >'
	     read(5,*)r_offset	     
	   write(6,'(a)')' Track electrons backwards? (y/n) >'
	   read(5,'(a)')ans
	   if(ans .eq. 'y' .or. ans .eq. 'Y') then
             backwards = .true.
	     write(6,'(a)')' Enter vertical incident angle >'
	     read(5,*)tht_in
!	     write(6,'(a)')' Enter radial offset >'
!	     read(5,*)r_offset	     
	     write(6,'(a)')' Enter incoming spin angles theta, phi (beam along x) >'
	     read(5,*) theta_s,phi_s   
	   end if
       	endif
	write(6,'(a)')' View along x,y,z or oblique? (x/y/z/o) >'
	read(5,'(a)')ans
	if(ans .eq. 'x' .or. ans .eq. 'X')view_x = .true.
	if(ans .eq. 'y' .or. ans .eq. 'Y')view_y = .true.
	if(ans .eq. 'z' .or. ans .eq. 'Z')view_z = .true.
c
c	Use Oxford field map or do uniform field test case
c
	if(.not. uniform_test)then	!read in numerical field map
	Bscale = B_max/5.1d0
	lun = 15
	   open(unit=lun,name='field_map.dat',status='old')
!C	   open(unit=lun,name='FIELD_MAP.DAT',status='old')
	   do ir=1,51
	     rr(ir) = 2.d0*float(ir-1)
	     zz(ir) = 2.d0*float(ir-1)
	     do iz=1,51
	       read(lun,*)xx,xx,B_field_z(iz,ir),B_field_r(iz,ir),xx,xx,xx
	       B_field_z(iz,ir) = B_field_z(iz,ir)*Bscale
	       B_field_r(iz,ir) = B_field_r(iz,ir)*Bscale
	     enddo
	   enddo
	   close(unit=lun)
	else
	   do ir=1,51			!uniform 5T field over 26 cm in z
	     rr(ir) = 2.d0*float(ir-1)	! and 16 cm in r
	     zz(ir) = 2.d0*float(ir-1)
	     do iz=1,51
	        B_field_r(iz,ir) = 0.d0
	        if(rr(ir).le.16. .and. zz(ir).le.26.)then
		    B_field_z(iz,ir) = 5.0
	        else
		    B_field_z(iz,ir) = 0.d0
	        endif
	     enddo
	   enddo
	endif
c
c       calculate kinematic quantities
c
	e0_e = e0_e*1000.d0					!to MeV
!	ef_e  = e0_e/(1.+2.*e0_e*sind(tht_e/2.)**2/tmass) - ebind
	ef_e  = (e0_e - ebind/(1.0d0 - ebind/tmass))/
	1	(1.0d0+2.0d0*e0_e*sind(tht_e/2.0d0)**2/
	1	(tmass - ebind))
	if(ansq .eq. 'q' .or. ansq .eq. 'Q')then
	 W = pmass
	else
	W = W*1000.d0						!to MeV
	 K = (pmass**2 - W**2)/(2.0d0*pmass)
	 ef_e  = (e0_e + K)/
	1	(1.0d0+2.0d0*e0_e*sind(tht_e/2.0d0)**2/pmass)
	end if
!	qvec  = sqrt((e0_e-ef_e)*2.*tmass + (e0_e-ef_e)**2)
	eloss =  e0_e - ef_e
	qsq = 4.0d0*e0_e*ef_e*sind(tht_e/2.0d0)**2
	if(ansp .eq. 'e' .or. ansp .eq. 'E')then
!	                               ! Moeller, E' can be about = m_e
	 tmass2 = tmass**2
	 qsq = -2.0d0*(tmass2 - e0_e*ef_e +
	1	sqrt(e0_e**2 - tmass2)*sqrt(ef_e**2 - tmass2)
	1	*cosd(tht_e))
	end if
	qvec  = sqrt(eloss**2+qsq)
	if(qvec.eq.0.d0) qvec=1.0d-08
	tht_q0 = ef_e*sind(tht_e)/qvec
	tht_q = -asind(tht_q0)
!	tht_q = -asin(tht_q)/pi*180.d0
	if(ansp.eq.'e') then ! Moeller
	 tht_q0 = sqrt(ef_e**2 - tmass2)*sind(tht_e)/qvec
	 tht_q = -asind(tht_q0)
!	 tht_q = -asin(sqrt(ef_e**2 - tmass2)*sind(tht_e)/qvec)*180.d0/pi
	end if	
	phi_max = phi_step*float(n_events)	
!	jphistep =  int(phi_step) +1	
c
c	set up initial "velocity" vector, 1st event is in (x,y) plane
c
	if(recoils)then			!tracking recoils (protons or deuterons)
	  e_tot = sqrt(qvec**2 + tmass**2)
	  beta     = qvec/e_tot
	  v_i(1) = beta*c*cosd(tht_q)
	  v_i(2) = beta*c*sind(tht_q)
!	  v_i(3) = beta*c*sind(tht_q)
	  v_i(3) = 0.d0
	  factor = 90.d0/e_tot			!gives ns^-1 from B in Tesla
	if (ansp.eq.'e') factor = -90.d0/e_tot	!gives ns^-1 from B in Tesla
	  gamma = e_tot/tmass
	elseif(backwards)then			!tracking incident beam
	  tgt_r = 0.0d0
	  beta = 1.0d0
	  e_tot = e0_e
C	  v_i(1) = -c
C	  v_i(2) = 0.d0
C	  v_i(3) = 0.d0
	  v_i(1) = -c*cosd(tht_in)
!	  v_i(2) = -c*sind(tht_in)
	  v_i(2) = 0.d0
	  v_i(3) = -c*sind(tht_in)
!	  factor = -90.d0/e_tot
	  factor = -90.d0/e_tot*charge
	  gamma = e_tot/emass
	  beta = sqrt(gamma**2 - 1.0d0)/gamma
	phi_step = 0.d0	! no azimuthal angle for incident beam
	phi_max = 0.d0
	else					!tracing scattered e-
	  e_tot = ef_e
!	  e_tot = e0_e
!	  beta = 1.0d0
	  v_i(1) = c*cosd(tht_e)
	  v_i(2) = c*sind(tht_e)
!	  v_i(3) = c*sind(tht_e)	! restore phi_lab dependence; must be zero for r-offset
	  v_i(3) = 0.d0
	if(r_offset.ne.0.d0) v_i(3) = 0.d0 ! no out-of-plane velocity
	  factor = -90.0d0/e_tot*charge
	  gamma = e_tot/emass
	  beta = sqrt(gamma**2 - 1.0d0)/gamma
!	  if(abs(phi_b-tht_e).gt.90.0d0) back_comp = .true.
	  if(abs(phi_b-tht_e).gt.90.0d0.and.
	1	abs(phi_b-tht_e).lt.270.0d0) back_comp = .true.
	endif
c
c	set up plot limits
c
	if(.not. backwards)then
	  xmin = -1.2d0*tgt_r
	  xmax = 1.5d0*v_i(1)*t_max
	  ymin = -1.2d0*v_i(2)*t_max-2.0d0
	  ymax = 1.2d0*v_i(2)*t_max+2.0
!	  zmin = -1.2d0*v_i(2)*t_max*sind(phi_max)- 4.0d0
!	  zmax = 1.2d0*v_i(2)*t_max*sind(phi_max)+ 4.0d0
!	  zmin = -1.2d0*v_i(3)*t_max- 4.0d0
!	  zmax = 1.2d0*v_i(3)*t_max+ 4.0d0
	  zmin = -100.d0
	  zmax = 100.d0
	else
	  xmax = 1.2d0*tgt_r
	  xmin = 1.5d0*v_i(1)*t_max
!	  ymin = -40.d0
!	  ymax = 40.d0
!	  zmin = -10.d0
!	  zmax = 10.d0
	  ymin = -100.d0
	  ymax = 100.d0
	  zmin = -100.d0
	  zmax = 100.d0
	endif
	
	open(unit=10,name='bdl.dat',status='unknown')
!	open(unit=11,name='bdl.sav',status='unknown')
	open(unit=7,name='path_map.top',status='unknown')
	open(unit=8,name='brz.dat',status='unknown')
	call header(xmin,xmax,ymin,ymax,zmin,zmax)
c
	t_step = t_max/float(n_steps)		!time step for solving ODE's
	n_track = 0
c

!	do j=1,2				!for left/right scatterings
!	jphimax = int(phi_max/jphistep) +1
!	 do jphi_var=-jphimax,jphimax,jphistep	!azimuthal angle loop
!	 do jphi=-n_events,n_events		!azimuthal angle loop

	 do jphi=1,n_events			!azimuthal angle loop

	 prcs(1) = sind(theta_s)*cosd(phi_s)
	 prcs(2) = sind(theta_s)*sind(phi_s)
	 prcs(3) = cosd(theta_s)
	  n_track = n_track + 1
!	  if(j .eq. 1)then
!	  if((n_track - 2*(n_track/2)).eq. 1)then
	    phi = jphi*phi_step
!	    phi = jphi_var*phi_step
!	  else
!	    phi = 180.d0 + jphi_var*phi_step
!	  endif
c
c	  set up initial conditions, including actual velocity vector
c
	  t = 0.d0
	  y(1) = origin*(2.*tgt_r*ran2(iseed) - tgt_r)	!x
	  y(2) = v_i(1)				!dx/dt
!	  y(3) = 0.d0				!y
	  y(3) = r_offset*cosd(phi)		!y
	  y(4) = v_i(2)*cosd(phi)		!dy/dt
!	  y(4) = v_i(2)			!dy/dt
!	  y(5) = 0.d0				!z
	  y(5) = r_offset*sind(phi)		!z
!  	  y(6) = v_i(3)			!dz/dt
  	  y(6) = v_i(3)*sind(phi)		!dz/dt
!c	  y(6) = v_i(2)*sind(phi)		!dz/dt
	  y(7) = 0.d0				!integral B_x-dl
	  y(8) = 0.d0				!integral B_y-dl
	  y(9) = 0.d0				!integral B_z-dl
c
	  do i=1,3				!save initial velocity
	   v_i_save(i) = y(2*i)
	  enddo
c
	  write(7,*)y(1),y(3),y(5)
cxxx	  write(8,*)y(2),y(4),y(6)		!to write out velocities
c
 	  do ii=1,n_steps			!step thru time
	    call derivs(t,y,dydx)
      do 11 i=1,9
        ywatch(i)=y(i)
        dywatch(i)=dydx(i)
11    continue
	    call rk4(y,dydx,9,t,t_step,yout)	!solve diff. eqs.
	    do i=1,9
		y(i) = yout(i)
	    enddo
	    write(7,*)yout(1),yout(3),yout(5)
cxxx	    write(8,*)yout(2),yout(4),yout(6)	!to write out velocities
	    t = t + t_step
c
c	intermediate Bdl for precession
c
	  theta_p = 171.887d0*yout(7)/e_tot/beta	! Bdl due to Bz
	  phi_p = 171.887d0*
	1	sqrt(yout(8)**2+yout(9)**2)/e_tot/beta! Bdl due to Br
	 
c
c	rotate spin to new precessed direction
c
	  call rot(prco,prcs,theta_p,phi_p)

	  do ip = 1,3
	    prcs(ip) = prco(ip)
	  enddo

	  enddo

	  write(7,*)'join solid'

	B_axis(1) = 30.d0*t_max
	B_axis(2) = 0.d0
	B_axis(3) = 0.d0
	call rot(B_lab,B_axis,0.d0,-phi_b)
	write(7,*)0., 0., 0.
	write(7,*)B_lab(1), B_lab(2), B_lab(3)
	  write(7,*)'join dots'

	BETA_axis(1) = 400.d0
	BETA_axis(2) = 0.d0
	BETA_axis(3) = 0.d0
	call rot(BETA_lab,BETA_axis,0.d0,-tht_e)
	write(7,*)0., 0., 0.
	write(7,*)BETA_lab(1), BETA_lab(2), BETA_lab(3)
	  write(7,*)'join dotdash'
c
c	  calculate total deflection angle based on integral B-dl
c
	  bdl = 0.d0
	  do k=7,9
	   bdl = bdl + yout(k)**2		!B_x^2+B_y^2+B_z^2
	  enddo
	  bdl = sqrt(bdl)
	  theta_bdl = 171.887d0*bdl/e_tot/beta	!deflection angle
	  theta_pre = gamma*g22*theta_bdl		!precession   "
c
c	  also calculate it from v_i-dot-v_f
c
	  arg1 = v_i_save(1)*yout(2) + v_i_save(2)*yout(4) 
     #	       + v_i_save(3)*yout(6)
	  arg2 = v_i_save(1)**2 + v_i_save(2)**2 + v_i_save(3)**2
	  arg3 = yout(2)**2 + yout(4)**2 + yout(6)**2
	  theta_if = acosd(arg1/sqrt(arg2*arg3))
c
	  if(theta_if .ne. 0.d0)ratio = theta_bdl/theta_if
	forb = 'f'
	if(backwards) forb = 'b'
	  write(10,1000)n_track,bdl,theta_if,theta_bdl,ratio
	q_sign = '-'
	if(charge.lt.0.d0) q_sign = '+' ! there is an extra (-) with factor: default is e-
	  write(10,1001)e0_e,tht_e,forb,
	1	tht_b,phi_b,anse,q_sign,ansp,
	1	tht_q, phi,B_max,W,e_tot,t_max,
	1	prcs(1),prcs(2),prcs(3)
!	1	theta_pre, phi,e_tot
1001	format(/,' Beam E ',g13.7,'theta ',g13.7,' backw.? ',a3,
	1	/,' B field theta, phi ',2g13.7,'e/recoil ',a2,a1,
	1	'; p/d ',a3,/,' theta_q ',g13.7,
!	1	'p/d ',a3,/,' theta_precession ',g13.7,
	1	'starting phi ',g13.7,'  B[T] ',g13.7,/,
	1	' W ',g13.7,' Final E',g13.7,'  t_max[ns] ',g13.7,/,
	1	' precession ',3g13.7)
c
c	Donal'sc	D
c
   	theta_bend = atand(sqrt(yout(2)**2 + yout(4)**2)/yout(6))
 	phi_bend = atand(yout(4)/yout(2))
	if(recoils) then
 	phi_bend = phi_bend - abs(tht_q) - abs(tht_e)
 	end if

c   
	if(backwards) then
   	write(10,1002)90-theta_bend,phi_bend
	else
   	write(10,1002)90+theta_bend,phi_bend-tht_e
!   	write(10,1002)(90+theta_bend)*charge,phi_bend-tht_e
	end if
1002	format(/,' Theta_bend = ',g13.7,' Phi_bend (w.r.t. theta) = ',g13.7)
!   	write(10,'(/,a)')' (Multiply by -1 for electrons)'
   
	 enddo					!end of 'phi' loop
!	 enddo					!end of 'phi_var' loop
!	enddo					!end of 'j=1,2' loop
c

!	call boundary(tgt_r,field_r)		!draw circle
1000	format(/,t10,'Track #: ',i3,t40,'Integral B_dl = ',f8.3,
     #	       ' T-cm',/,t10,'angle(v_i,v_f) = ',f9.3,
     #	       t40,'angle_Bdl = ',f9.3,' degrees',/,
     #	       t10,'Bdl/v_if  = ',f9.4)
	stop
	end	
c
c-------------------------------------------------------------------------------
c	The derivatives of (x,dx/dt,y,dy/dt,z,dz/dt) wrt. time
c
	subroutine derivs(x,y,dydx)
	implicit real*8(a-h,o-z)
	dimension y(9),dydx(9),w_f(3)
	common/conversion_factor/factor
c
	call cyclotron_freq(y,w_f)
c
c	These are just the velocities
c
	dydx(1) = y(2)
	dydx(3) = y(4)
	dydx(5) = y(6)
!	write(11,100) dydx(1),dydx(3),dydx(5)
100	format(//4g20.13)
c
c	This is just v_vec X omega_vec
c
	dydx(2) = y(4)*w_f(3) - y(6)*w_f(2)
	dydx(4) = y(6)*w_f(1) - y(2)*w_f(3)
	dydx(6) = y(2)*w_f(2) - y(4)*w_f(1)
!	write(11,101) dydx(2),dydx(4),dydx(6),w_f(1),w_f(2),w_f(3)
101	format(3g20.13)
c
c	These are the components of B_perp
c
	dydx(7) = dydx(2)/factor
	dydx(8) = dydx(4)/factor
	dydx(9) = dydx(6)/factor
!	write(11,100) dydx(7),dydx(8),dydx(9),factor

c
	return
	end
c------------------------------------------------------------------------------ 
	subroutine cyclotron_freq(yp,w)
c
c	Calculates v_vec X B_vec for arbitrarily orientated field
c
	implicit real*8(a-h,o-z)
!	real*4 tht
	logical*1 backwards,back_comp
	dimension B_field_z(51,51),B_field_r(51,51),zz(51),rr(51)
	dimension w(3),yp(9),coord_lab(3),coord_prime(3),w_p(3)
c
	common/field_angles/tht_b,phi_b
	common/field_strength/B_field_z,B_field_r,zz,rr
	common/conversion_factor/factor
        common/back/backwards,back_comp	! "backward" component in v_i,
!						  because phi_b+ tht_e > 90 degrees.	        	

!	parameter (save_rc=-1.d0) ! typo
	data save_r_c /-1.d0/

	external sind,cosd,tand,asind,acosd,atand
c
	coord_lab(1) = yp(1)		!(x,y,z)_LAB
	coord_lab(2) = yp(3)
	coord_lab(3) = yp(5)
c
c	rotate to coordinates with z' along field direction
c
	call rot(coord_prime,coord_lab,tht_b,phi_b)

c
c	compute azimuthal angle, go thru a bunch of stuff to ensure
c	proper quadrant
c
	isigns = 1
	isignc = 1
  	if(coord_prime(1).ne. 0.d0)tht = 
     #		dabs(atand(coord_prime(2)/coord_prime(1)))
!  	if(coord_prime(1).ne. 0.d0)tht = 
!     #		dabs(atan(coord_prime(2)/coord_prime(1))*180/pi)
	if(coord_prime(1) .lt. 0.d0)then
!	   tht = tht + 180.d0
	   isignc = -1
	   if(coord_prime(2) .lt. 0.d0)then
		isigns = -1
	   end if
  	elseif(coord_prime(2).lt.0.d0 .and. coord_prime(1).gt. 0.d0)then
!	   tht = tht + 360.d0
	   isigns = -1
  	elseif(coord_prime(2).lt.0.d0 .and. coord_prime(1).eq. 0.d0)then
	   tht = 270.d0
  	elseif(coord_prime(2).gt.0.d0 .and. coord_prime(1).eq. 0.d0)then
	   tht = 90.d0
  	elseif(coord_prime(2).eq.0.d0 .and. coord_prime(1).eq. 0.d0)then
	   tht = 0.d0
	endif

	z_c = abs(coord_prime(3))
	r_c = sqrt(coord_prime(1)**2 + coord_prime(2)**2)


c
c	compute field in target coordinate system
c
	b_z = rinterp2d(zz,rr,B_field_z,z_c,r_c,51,51,2.0d0,2.0d0)
	b_r = rinterp2d(zz,rr,B_field_r,z_c,r_c,51,51,2.0d0,2.0d0)

c	if(phi_b.gt.90.d0.and.phi_b.lt.270.d0) b_r = -b_r
c      This was old way fo reversing the field, does not always work

! *****  reverse field direction for angles < 270 deg and > 90 deg. *****
! *****  for tracking the beam  *****

	if(backwards) then
	if(phi_b.gt.90.d0.and.phi_b.lt.270.d0) b_r = -b_r
	isignc = -isignc
	isigns = -isigns
	end if

! *****  for tracking recoil (scattered) particles  *****
! *****  back_comp is initialized as .false. in main, set .true. if sum  *****
! *****  of angles abs(phi_b-tht_e) is > 90 and < 270, also in main initialization  ***** 

	if(back_comp) then
	isignc = -isignc
	isigns = -isigns
	end if

!	if(r_c.gt.0.d0) then
!	w_p(1) = factor*b_r*coord_prime(1)/r_c		!local cyclotron freq.
!	w_p(2) = factor*b_r*coord_prime(2)/r_c
!!	write(6,*)w_p(1), w_p(2), coord_prime(1)/r_c
!	else
!	w_p(1) = 0.d0
!	w_p(2) = 0.d0
!	end if

	w_p(1) = factor*b_r*cosd(tht)*isignc		!local cyclotron freq.
	w_p(2) = factor*b_r*sind(tht)*isigns
	w_p(3) = factor*b_z

!	w_p(3) = factor*b_z
c
c	transform back to LAB
c
	call roti(w,w_p,tht_b,phi_b)
c
!	write(7,*)r_c,z_c,b_r,b_z,w

	if(save_r_c.ne.r_c) then	! to write only once per deriv call.
	write(6,*)r_c,z_c,b_r,b_z
	write(8,*)r_c,z_c,b_r,b_z
	save_r_c = r_c
	end if

	return
	end
C------------------------------------------------------------------------------ 
C     2-Dimensional Linear Interpolation Routine:                               
C                                                                               
C       Calculates F(X0,Y0) given array F(X,Y) by two successive                
C       interpolations, first in X and then in Y.                               
C                                                                               
C       Assumes uniform spacing of array in X and Y.                            
C------------------------------------------------------------------------------ 
c                                                                               
      real*8 function rinterp2d(x,y,f,x0,y0,nx,ny,delx,dely)          
	implicit real*8(a-h,o-z)
      dimension x(51),y(51),f(51,51)                                        
      logical*1 uniform_test
      common/field_type/uniform_test
c                                                                               
      if(uniform_test)then			!outside uniform region?
	if(x0 .gt. 26. .or. y0 .gt. 16.)then
	   rinterp2d = 0.d0
	   return
	endif
      endif
c
      i = int((x0-x(1))/delx) + 1                                              
      j = int((y0-y(1))/dely) + 1                                              
      if((i+1.gt.nx).or.(i.lt.1).or.(j+1.gt.ny).or.(j.lt.1))then                
        rinterp2d = 0.d0
        return                                                                  
      endif                                                                     
c                                                                               
      rintx1 = ((x0-x(i))/(x(i+1)-x(i)))*(f(i+1,j)-f(i,j))                      
     #                  + f(i,j)                                                
      rintx2 = ((x0-x(i))/(x(i+1)-x(i)))*(f(i+1,j+1)-f(i,j+1))                  
     #                  + f(i,j+1)                                              
      rinterp2d = ((y0-y(j))/(y(j+1)-y(j)))*(rintx2-rintx1) + rintx1            
c                                                                               
      return                                                                    
      end                                                                       
c----------------------------------------------------------------------------
c	Fourth-order Runge-Kutta from Numerical Recipes book
c
      subroutine rk4(y,dydx,n,x,h,yout)
	implicit real*8(a-h,o-z)
!      parameter (nmax=10)
!      dimension y(9),dydx(9),yout(9),yt(nmax),dyt(nmax),dym(nmax)
      dimension y(9),dydx(9),yout(9),yt(10),dyt(10),dym(10)

      hh=h*0.5d0
      h6=h/6.d0
      xh=x+hh
      do 11 i=1,n
        yt(i)=y(i)+hh*dydx(i)
11    continue
!	write(11,100) dyt
      call derivs(xh,yt,dyt)
!	write(11,100) dyt
      do 12 i=1,n
        yt(i)=y(i)+hh*dyt(i)
12    continue
!	write(11,101)
101	format(' Call deriv')
102	format(' Last deriv')
! 	write(11,100) dym
      call derivs(xh,yt,dym)
! 	write(11,100) dym
      do 13 i=1,n
        yt(i)=y(i)+h*dym(i)
        dym(i)=dyt(i)+dym(i)
13    continue
!	write(11,101)
!	write(11,100) dyt
! 	write(11,100) dym
100	format(4g20.13)

      call derivs(x+h,yt,dyt)
!	write(11,102)
!	write(11,100) dyt
      do 14 i=1,n
        yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.d0*dym(i))
14    continue
      return
      end
c------------------------------------------------------------------------
	subroutine boundary(tgt_r,field_r)
	implicit real*8(a-h,o-z)
c
c	draw circles corresponding to target material and field boundaries
c
	write(7,*)'set polar'

	do tht=0,360,10				!field
	   write(7,*)field_r,tht
	enddo
	write(7,*)'join dots'

	do tht=0,360,10				!target
	   write(7,*)tgt_r,tht
	enddo
	write(7,*)'join dotdash'

	return
	end
c------------------------------------------------------------------------------
c	random number routine from Numerical Recipes book
c	gives a uniform deviate on [0,1) with no sequential correlations
c
	real*8 function ran2(idum)
	implicit real*8(a-h,o-z)
	dimension ir(97)
	data iff/0/
	m=714025
	ia=1366
	ic=150889
	rm=1./m
	if(idum .lt. 0 .or. iff .eq. 0)then
	    iff = 1
	    idum = mod(ic-idum,m)
	    do j=1,97
		idum = mod(ia*idum+ic,m)
		ir(j) = idum
	    enddo
	    idum = mod(ia*idum+ic,m)
	    iy = idum
	endif
	j = 1 + (97*iy)/m
	if(j .gt. 97 .or. j .lt. 1)pause
	iy = ir(j)
	ran2 = iy*rm
	idum = mod(ia*idum+ic,m)
	ir(j) = idum
	return
	end
c
c-----------------------------------------------------------------------------
c	header stuff for plot file
c
	subroutine header(xmin,xmax,ymin,ymax,zmin,zmax)
	implicit real*8(a-h,o-z)
	logical*1 view_x,view_y,view_z
	common/view/view_x,view_y,view_z
c
!	write(7,*)'set outline off'
	write(7,*)'set outline on'
	write(7,*)'set labels size 4.5'
	write(7,*)'set order x y z'
	write(7,*)'set three scrd 14'
	if(view_x)then
	  write(7,*)'set three theta 0.'
	  write(7,*)'set three phi 90.'
	elseif(view_y)then
	  write(7,*)'set three theta 90.'
	  write(7,*)'set three phi 90.'
	elseif(view_z)then
	  write(7,*)'set three theta 0.'
	  write(7,*)'set three phi 0.'
	else
	  write(7,*)'set three theta 20.'
	  write(7,*)'set three phi 100.'
	endif
	write(7,*)'set limits x',xmin,xmax
	write(7,*)'set limits y',ymin,ymax
	write(7,*)'set limits z',zmin,zmax
!	write(7,*)1.5*xmin, 0., 0.
!	write(7,*)1.5*xmax, 0., 0.
	write(7,*)1.1*xmin, 0., 0.
	write(7,*)1.2*xmax, 0., 0.
	write(7,*)'set pattern 0.05 0.05'
	write(7,*)'join pattern'
	write(7,*)'plot axes'
c
	return
	end
c
c------------------------------------------------------------------------------ 
	subroutine rot(u,v,t,p)
c
c	input vector v is rotated to coordinate system whose z axis is
c	specified by spherical coordinates (t,p), i.e., (theta,phi)
c
	implicit real*8(a-h,o-z)
	dimension u(3),v(3)
	external sind,cosd,tand,asind,acosd,atand
c
	ct=cosd(t)
	st=sind(t)
	cp=cosd(p)
	sp=sind(p)
c
	u(1) = cp*ct*v(1) + ct*sp*v(2) - st*v(3)
	u(2) = -sp*v(1) + cp*v(2)
	u(3) = cp*st*v(1) + st*sp*v(2) + ct*v(3)
c
	return
	end

c------------------------------------------------------------------------------ 
	subroutine roti(u,v,t,p)
c
c	Inverse of rotation performed by ROT
c
	implicit real*8(a-h,o-z)
	dimension u(3),v(3)
	external sind,cosd,tand,asind,acosd,atand
c
	ct=cosd(t)
	st=sind(t)
	cp=cosd(p)
	sp=sind(p)
c
	u(1) = cp*ct*v(1) - sp*v(2) + cp*st*v(3)
	u(2) = ct*sp*v(1) + cp*v(2) + st*sp*v(3)
	u(3) = -st*v(1) + ct*v(3)
c
	return
	end
c------------------------------------------------------------------------------ 
	real*8 function sind(x)
	real*8 pi,x,x0
	pi=3.141592654d0
	x0 = pi*x/180.0d0
	sind = sin(x0)
	return
	end
c------------------------------------------------------------------------------ 
	real*8 function cosd(x)
	real*8 pi,x,x0
	pi=3.141592654d0
	x0 = pi*x/180.0d0	
	cosd = cos(x0)
	return      
	end
c------------------------------------------------------------------------------ 
	real*8 function tand(x)
	real*8 pi,x,x0
	pi=3.141592654d0
	x0 = pi*x/180.0d0
	tand = tan(x0)
	return
	end
c------------------------------------------------------------------------------ 
	real*8 function asind(x)
	real*8 pi,x,x0
	pi=3.141592654d0
	x0 = x
	asind = asin(x0)*180.d0/pi
	return
	end
c------------------------------------------------------------------------------ 
	real*8 function acosd(x)
	real*8 pi,x,x0
	pi=3.141592654d0
	x0 = x
	acosd = acos(x)*180.d0/pi
	return
	end
c------------------------------------------------------------------------------ 
	real*8 function atand(x)
	real*8 pi,x,x0
	pi=3.141592654d0
	x0 = x
	atand = atan(x0)*180.d0/pi
	return
	end