C These routines, for magnetic field line tracing, where taken from C a field-line tracer developed for the Toffoletto-Hill '93 model. c ---------- predraw ---------------------------------------------- c this subroutine sets up the parameters for the program c er = field line tracer error, typically=1.0e-03 subroutine predraw(x1,x2,y1,y2,z1,z2,imax1,err1,err2) implicit double precision(a-h,o-z) real x1,x2,y1,y2,z1,z2,err1,err2 integer imax1 parameter (pi=3.14159265358) common /calrange/ xmin,xmax,ymin,ymax,zmax,zmin,imax common /errrange/ er1,er2 c ... for the calculation region ................................. C write(6,*) 'Enter the tracing parameters : ' C write(6,*) ' enter xmax, xmin, ymax, ymin, zmax, zmin ' C read(5,*) xmax,xmin,ymax,ymin,zmax,zmin C write(6,*) ' xmax = ',sngl(xmax),' xmin =',sngl(xmin) C write(6,*) ' ymax = ',sngl(ymax),' ymin =',sngl(ymin) C write(6,*) ' zmax = ',sngl(zmax),' zmin =',sngl(zmin) xmin=x1 xmax=x2 ymin=y1 ymax=y2 zmin=z1 zmax=z2 imax=imax1 C 10 write(6,*) ' enter the maximum tracing steps (imax<2000) :' C write(6,*) ' ( imax=0 : no limit on the number of steps) ' C read(5,*) imax C if (imax.gt.2000) goto 10 C write(6,*) ' imax =',imax C write(6,*) ' enter the tracing errors er1 and er2 ' C read(5,*) er1,er2 C write(6,*) ' er1 =',sngl(er1),' er2 =',sngl(er2) er1=err1 er2=err2 return end c ---------- end of predraw --------------------------------------- c ------------- drawline ----------------------------------------- subroutine drawline(xx0,yy0,zz0,dir,xx,yy,zz,idim,iimax) implicit double precision(a-h,o-z) real xx(idim),yy(idim),zz(idim) common /calrange/ xmin,xmax,ymin,ymax,zmax,zmin,imax common /errrange/ er1,er2 common /mpause/radiusmp,dnosemp,bexmp common /bstep/ dx1,dy1,dz1,bb,h,ep,er data bmin0 /0.0001/, grr /1.0/, h /0.25/,eps /0.01/ poten = 999.0 flen = 0.0 dt = 0.0 er = er1 ep = dir yymax=max(abs(ymax),abs(ymin)) zzmax=max(abs(zmax),abs(zmin)) rmax=sqrt(yymax*yymax+zzmax*zzmax) if(rmax.lt.20.0) rmax=20.0 dx1 = 0.0 dxold = 0.0 iimax = 1 iopen = 0 ineutr = 0 xm = xx0 ym = yy0 zm = zz0 xnmax = xx0 ynmax = yy0 znmax = zz0 xr = 0.0 yr = 0.0 zr = 0.0 xmp = 0.0 ymp = 0.0 zmp = 0.0 c ...... reset the array .............. do 900 i=1,idim xx(i) = 0.0 yy(i) = 0.0 zz(i) = 0.0 900 continue xx(1) = xx0 yy(1) = yy0 zz(1) = zz0 i = 2 c ........... field line loop ................................ 405 continue call trac(xm,ym,zm,xn,yn,zn) c field line length calculation ........................ dr = sqrt((xm-xn)**2+(ym-yn)**2+(zm-zn)**2) flen = flen + dr c flux tube volume calculation ......................... if(bb.ge.bmin0) dt = dt+dr/bb rrr = sqrt(xn*xn + yn*yn + zn*zn) rrr_old = sqrt(xm*xm + ym*ym + zm*zm) rr = sqrt(yn*yn+zn*zn) if(xn+bexmp.le.0.0)then rmod = rr else rmod = sqrt((xn+bexmp)**2+yn*yn+zn*zn) end if if(xm+bexmp.le.0.0)then rmod_old = sqrt(ym*ym+zm*zm) else rmod_old = sqrt((xm+bexmp)**2+ym*ym+zm*zm) end if c ........determine the farthest distance away from earth .... if(rrr.gt.rrr_old) then xnmax=xn ynmax=yn znmax=zn end if c ........ case: inside the earth ............................ if(rrr.lt.grr)then C call transsph(xn,yn,zn,gl,gp) C write(6,806) sngl(gl),sngl(gp) xx(i) = xn yy(i) = yn zz(i) = zn iopen = -1 rrr = sqrt(xn**2+yn**2+zn**2) if(rrr.eq.0.0) rrr=1.0 goto 404 end if c ...... case: crossed the magnetopause ....................... C if((rmod-radiusmp)*(rmod_old-radiusmp).lt.0.0)then C iopen = 1 C xmp = 0.5*(xm+xn) C ymp = 0.5*(ym+yn) C zmp = 0.5*(zm+zn) C write(6,800) sngl(xmp),sngl(ymp),sngl(zmp) C end if c ........ case: cross the neutral sheet ......................... if(dx1*dxold.lt.0.0)then if(xm+bexmp.lt.0.0.or.abs(zn).lt.0.5)then ineutr = 1 xr = xn yr = yn zr = zn br = bb write(6,811)sngl(xr),sngl(yr),sngl(zr) end if end if dxold = dx1 c ........ case: cross the equatorial plane ................... if(zm*zn.lt.0.)then write(6,801) sngl(xn),sngl(yn),sngl(zn) c if(iplot.eq.3) then c xeq=(xm+xn)/2.0 c yeq=(ym+yn)/2.0 c zeq=(zm+zn)/2.0 c endif end if c ........ case: exceeded the maximum step ......................... if((i.gt.imax).and.(imax.gt.0)) then write(6,*) ' exceeded max #steps before crossing the boundary' goto 404 endif c ........ case: exceeded limits ............................ if((xn.le.xmin).or.(yn.le.ymin).or. 1 (zn.le.zmin).or.(zn.ge.zmax).or. 2 (xn.ge.xmax).or.(yn.ge.ymax).or. 2 (rr.gt.rmax)) 3 then xx(i) = xn yy(i) = yn zz(i) = zn go to 404 end if c ........ case: field is too small ............................ if(bb.lt.bmin0) then write(6,*)' !! |B| << ',sngl(bmin0) goto 404 end if c ..... assign array values ....................... xx(i) = xn yy(i) = yn zz(i) = zn xm = xn ym = yn zm = zn i = i + 1 if(h.gt.2.0) h=2.0 if(abs(rmod-radiusmp).gt.0.1*radiusmp)er=er2 goto 405 c ........... end of field line loop ........................... 404 continue iimax = i - 1 c ....... begin writing out the results ..................... write(6,802) sngl(xn),sngl(yn),sngl(zn) write(6,803) sngl(xnmax),sngl(ynmax),sngl(znmax) write(6,804) iimax write(6,805) sngl(flen) write(6,809) sngl(dt) c ......... open field line potential ..................... if(iopen.eq.1)then C call calpot(xmp,ymp,zmp,poten) C write(6,807) sngl(poten) end if c ......... closed field line potential ..................... if(ineutr.eq.1.and.iopen.eq.-1) iopen=-2 if(iopen.eq.-2)then rrr1=sqrt(xr**2+yr**2+zr**2) if(rrr1.ne.0.0)then c call calpotc(xr,yr,zr,poten) C write(6,808) sngl(poten) end if end if 800 format(2x,27Hcrossed the magnetopause at,1x,3f10.4) 801 format(2x,27Hcrossed equatorial plane at,1x,3f10.4) 802 format(2x,17Hfinal position at,1x,3f10.4) 803 format(2x,20Hfarthest position at,1x,3f10.4) 804 format(2x,14Htotal points =,I6) 805 format(2x,19Hfield line length =,f10.4,2x,2HRE) 809 format(2x,'flux tube volume =',f14.8,' RE/nT') 806 format(2x,28Hclosed on Earth at glat = ,f8.3,2x, * 9H glong = , f8.3) 807 format(2x,27Hopen field line potential =,f12.6,' KV') 808 format(2x,29Hclosed field line potential =,f12.6,' KV') 811 format(2x,24Hcrossed neutral sheet at,1x,3f10.4) 810 format(2f12.4,3x,f12.6) 999 format(3f12.4) 409 continue return end c ------------ end of drawline ----------------------------- c --------------v trac v-------------------------------------- subroutine trac(xa,ya,za,xn,yn,zn) implicit double precision (a-h,o-z) common /bstep/ dx1,dy1,dz1,bf1,h,ep,er call point93(xa,ya,za,dx1,dy1,dz1,bf1) if(bf1.eq.0.)return 1 continue h1 = h * ep h2 = 0.5 * h1 xn1 = xa + h1*dx1 yn1 = ya + h1*dy1 zn1 = za + h1*dz1 xn2 = xa + h2*dx1 yn2 = ya + h2*dy1 zn2 = za + h2*dz1 call point93(xn2,yn2,zn2,dx2,dy2,dz2,bf2) if(bf2.eq.0)return xn3 = xn2 + h2*dx2 yn3 = yn2 + h2*dy2 zn3 = zn2 + h2*dz2 delx = xn1 - xn3 dely = yn1 - yn3 delz = zn1 - zn3 rtes1 = delx*delx rtes2 = dely*dely rtes3 = delz*delz rtes = sqrt(rtes1 + rtes2 + rtes3) c rtes = rtes/sqrt(xa**2+ya**2+za**2) if(rtes-er) 4,4,2 4 continue xn = xn3 - delx yn = yn3 - dely zn = zn3 - delz c bb = bf2 if(rtes-(0.25*er)) 7,6,6 7 h = 2.0*h 6 go to 5 2 continue h = 0.5*h go to 1 5 continue return end c ------------v end of trac v----------------------------------- c ------------- transsph -------------------------------- subroutine transsph(x,y,z,glat,gphi) implicit double precision(a-h,o-z) parameter (pi=3.14159265358) r=sqrt(x**2+y**2+z**2) rho=sqrt(x**2+y**2) t = atan3(z,rho) glat = t/pi*180. gphi = 0.0 p = atan3(y,x) gphi = p/pi*180. return end c -------- end of transsph ------------------------------------ c *************************************************************** c c should be followed by Subroutine point93(x,y,z,bx,by,bz,b) c c ***************************************************************