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