c *********************************************************************** c c The Toffoletto and Hill 1993 open magnetosphere model c Variant for driving without console input c c Fortran 77 code c Version 1.20 c last modified 4/15/94 c variant 8/8/96 c c c changes: 9/92 original written by Frank Toffoletto and Tom Hill c 8/93 incorprate a flexible magnetopause nose distance. c 12/93 rewritten by Cheng Ding and Frank Toffoletto c 2/94 use the newer version of Voigt 81 model c 2/94 Bn on the magpause depends on cos(phi)**7 c 4/94 change sub "standoff"->"nosedist", "nose"->"dayn", etc. c 8/96 all inputs now set in call to one subroutine c c ********************************************************************** c --------------t point93 t--------------------------------------------- subroutine point93(xn,yn,zn,dx,dy,dz,bt,ptl,ptl0) implicit double precision(a-h,o-z) common /batch/ btail,z11,z21,ib common /mpause/radmp,dnosemp,bexmp common /block5/ csph,snph,qcs,qsn,csps,snps common /block3/ h0, alfa, n common /fimf/ himf, cdipol, cimf, xi common /rotate/ rot,fac common /imfrot/ rotate common /angdip/ tilt data tiny /1.0d-10/ data cdipol /0.0/, cimf /0.0/,bimf /0.0/, dir /0.0/ if(radmp.eq.0.0)then write(6,*)' magnetopause radius = 0.0 in point93 ' stop end if if(bexmp.eq.0.0)then write(6,*)' bexmp = 0.0 in point93' stop end if if(dnosemp.eq.0.0)then write(6,*)' magnetopause nose distance = 0.0 in point93' stop end if bx = 0.0 by = 0.0 bz = 0.0 c --- transform from GSM to model coordinates ------------- call coord(xn,yn,zn,bx,by,bz,-1) if(abs(xn-bexmp).lt.tiny) xn=bexmp if(abs(yn).lt.tiny) yn=0.0 if(abs(zn).lt.tiny) zn=0.0 rad = sqrt((xn-bexmp)**2 + yn*yn + zn*zn) if(xn-bexmp.gt.0.0) rad = sqrt(yn*yn+zn*zn) if(rad.lt.radmp)then c ............ internal field ................. iswit = 4 xgsm = -xn ygsm = -yn zgsm = zn call bvoigt(dnosemp,radmp,tilt,cdipol,cimf,bimf,dir, * iswit,xgsm,ygsm,zgsm,bxgsm,bygsm,bzgsm) bx = -bxgsm by = -bygsm bz = bzgsm c............. normal component field .................. if(fac.ne.0.)then iswit = 2 call norm(dnosemp,radmp,iswit,xn,yn,zn,bnx2,bny2,bnz2) call normp(dnosemp,radmp,xn-bexmp,yn,zn,bnx,bny,bnz,rat,iswit, 1 ptl,ptl0) bx = bx - fac*(bnx+bnx2) by = by - fac*(bny+bny2) bz = bz - fac*(bnz+bnz2) end if c ........... additional tail field ....................... if(btail.ne.0.)then iswit = 2 call tailst(dnosemp,radmp,iswit,xn,yn,zn,bx1,by1,bz1) if(rat.eq.0.0) rat=1.0 bx1 = - bx1*btail*rat by1 = - by1*btail*rat bz1 = - bz1*btail*rat bx = bx + bx1 by = by + by1 bz = bz + bz1 end if endif if(rad.ge.radmp)then c ............ magnetosheath field ................ yn1 = yn*cos(rotate) + zn*sin(rotate) zn1 = zn*cos(rotate) - yn*sin(rotate) call bint(xn-bexmp,yn1,zn1,bx,by1,bz1) by = by1*cos(rotate) - bz1*sin(rotate) bz = bz1*cos(rotate) + by1*sin(rotate) c ............ normal component field ............... if(fac.ne.0.)then call normp(dnosemp,radmp,xn-bexmp,yn,zn,bnx,bny,bnz,rat,2, 1 ptl,ptl0) bx = bx - fac*bnx by = by - fac*bny bz = bz - fac*bnz end if end if c --- transform to GSM coordinates ------------- call coord(xn,yn,zn,bx,by,bz,1) c ... the magnitude of the total B field .............. bt = sqrt(bx*bx+by*by+bz*bz) if(abs(bt).le.tiny)then write(6,*) ' bt = ',bt return end if dx = bx/bt dy = by/bt dz = bz/bt return end c ------------t end of point93 t--------------------------------- c ----------t setupth93 t---------------------------------------------- subroutine setupth93(dens,vsw,bxi,byi,bzi,tlt,ptl,ptl0,pz11,pz21) implicit double precision(a-h,o-z) common /imfb/ bimfx,bimfy,bimfz common /imfrot/ rotate common /mpause/ radius1,dnose1,bex1 common /angdip/ tilt common /batch/ btail,z11,z21,ib common /stretch/ str,rlam,tlam,y0 common /tailf/ bt,tw common /block1/ bex,r,q,pi,psi common /block3/ h0, alfa, n common /fimf/ himf, cdipol, cimf, xi common /rotate/ rot,fac data pi /3.14159265358/ data tilt/0.0/, radius/20.0/ c ... read in the IMF information ............................ c write(6,*) ' ' c write(6,*) 'Please enter the IMF information !' c write(6,*) ' enter solar wind speed (km/s) and density (part/cc) ' c read(5,*) vsw,dens c write(6,*) ' vsw =',sngl(vsw),' km/s, density =',sngl(dens),' /cc' c write(6,*) ' enter IMF Bx,By,Bz (nT) ' c read(5,*) bimfx,bimfy,bimfz bimfx=bxi bimfy=byi bimfz=bzi c write(6,*) ' bimfx=',sngl(bimfx),' bimfy=',sngl(bimfy), c * ' bimfz=',sngl(bimfz) bimf = sqrt(bimfx*bimfx+bimfy*bimfy+bimfz*bimfz) if(bimf.eq.0.0)then write(6,*) ' IMF B=0 gives a closed magnetosphere ! ' factor = 0.0 rot = 0.0 goto 101 end if c ... determine imf clock angle (bx is not used) ................ if(bimfz.eq.0.0.and.bimfy.eq.0.0.and.bimfx.ne.0.0)then c write(6,*) ' IMF by,bz cannot be zero while bx is nonzero !' c write(6,*) ' Sorry. Future versions may handle this !' stop end if if(bimfz.ne.0.0)then rot = atan2(-bimfy,-bimfz)/pi*180.0 else rot = -90.0 if(bimfy.lt.0.0) rot = 90.0 end if c write(6,*) ' IMF clock angle (from southward) =',sngl(rot) c ... factor for normal component and tail field ................. factor = cos(rot/2./180.*pi)**2 * bimf/6.0 *vsw/400.0 101 continue c ... calculate magnetopause nose distance .......................... call nosedist(vsw,dens,bimfy,bimfz,dnosemp) c write(6,*) ' magnetopause nose distance =',sngl(dnosemp),' Re' bex = radius-dnosemp bex1 = bex radius1 = radius dnose1 = dnosemp c write(6,*) ' ' c write(6,*) ' inputting magnetosheath data ......' call init c ... initialize the magnetosphere parameters ..................... c write(6,*) ' ' c write(6,*) 'Enter the information about inner magnetosphere : ' c write(6,*) ' enter the tilt angle of Earth dipole ' c write(6,*) ' ( + --> Northern hemisphere summer ) ' c read(5,*) tilt tilt=tlt c write(6,*) ' tilt =',sngl(tilt) c ... startup gehavo ................................................ alfa = 0.2 c write(6,*) ' calculating coeffs for stretching factor = ', c & sngl(alfa) call bvoigt(dnosemp,radius,tilt,cdipol,cimf,bimf,dir, * 1,xgsm,ygsm,zgsm,bxgsm,bygsm,bzgsm) call bvoigt(dnosemp,radius,tilt,cdipol,cimf,bimf,dir, * 2,xgsm,ygsm,zgsm,bxgsm,bygsm,bzgsm) c ... startup norm .......................................... c write(6,*) ' ' c used to be 1.5 here c fac = 1.5*factor fac=2.1875*factor rotate = rot/180.*3.14159 if(fac.ne.0.)then c write(6,*) 'Norm component ...' c write(6,*) ' average magnitude =',sngl(fac) rotml = rot/2.0 call norm(dnosemp,radius,1,rotml,0.0d0,0.0d0,dbx,dby,dbz) call normp(dnosemp,radius,rotml,0.0d0,0.0d0,dbx,dby,dbz,a,1, 1 ptl,ptl0) end if c ... startup tail ............................................. if(factor.ne.0.0) then c used to be 20.0 here btail = 30.0*factor c write(6,*) ' ' c write(6,'(a)') 'Tail field parameters ' c write(6,'(a)') ' enter plasma sheet thickness z11,z21 (eg.1,9)' tw = rot/2.0 stretch=0.0 rlam = 0.0 tlam = 0.0 y0 = 0.0 c read(5,*) z11,z21 z11=pz11 z21=pz21 c write(6,*)' tail field strength (nT) =',sngl(btail) c write(6,*)' plasma sheet thickness =',sngl(z11),sngl(z21) c write(6,*)' tail twist (= IMF clock angle/2) =',sngl(tw) ib=1 call tailst(dnosemp,radius,1,x,y,z,bx,by,bz) else btail=0.0 endif c ... input the potential .................................... C write(6,*) ' ' call potin(factor,rot/2.0) return end c --------------t end of setupth93 t------------------------------ c ------------t nosedist t---------------------------------------- c this subroutine gives the nose distance of the magnetopause c as a function of solar wind speed, density, and IMF direction c based of the Toffoletto and Voigt paper, AGU spring 1993 . c created 8/93 c c variables definition c c vsw - solar wind speed in km/s c den - solar wind density in particles/cc c ram - ram pressure in P c bimfx,bimfy, c bimfz,bimf - imf field in nt c kappa - Spreiter correction to ram pressure c f - compression factor of CF currents at magnetopause c alpha - interconnection field at nose for southward IMF c dm - dipole moment in nT c mp - mass of proton in kg c nt - a nanotesla c dnosemp - nose distance of magnetopause in RE (output) c alp - merging line angle c a,b,c,disc - variables used in calculation c br - ring current magnitude at nose (estimate) in nT c ---------------------------------------------------------- subroutine nosedist(vsw,den,bimfy,bimfz,dnosemp) implicit none real*8 vsw,den,ram,dnosemp real*8 bimfy,bimfz,bimf, br real*8 kappa,f,alpha, alp,a,b,c,disc real*8 dm,mp,nt,pi,mu0 parameter (pi=3.14159265358) data kappa /0.88/ data f /2.5/ data alpha /7.0/ data dm /3.1e-5/ data mp /1.67e-27/ data nt /1.0e-9/ data br /13e-9/ mu0 = 4.0*pi*1.0e-7 ram = kappa*den*1.0e6*mp*(vsw*1e3)**2 C write(6,*)' SW ram pressure =',sngl(ram)*1e09,' nP' ram = ram*2.0*mu0 bimf = sqrt(bimfy**2+bimfz**2) if(bimf.eq.0.0)then alp = 0.0 else alp = (1.0-bimfz/bimf)/2.0 end if a = ((alpha*alp)**2-1.0)*(bimf*nt)**2 * -ram+br**2 * -2.0*alpha*alp**1.5*br*bimf*nt b = -2.0*f*dm*bimf*nt*alp**1.5*alpha * +2.0*f*dm*br c = (f*dm)**2 disc = b**2-4.0*a*c if(disc.lt.0.)then write(6,*)' disc riminant in "nose" lt 0.0, aborting' stop end if dnosemp = -b-sqrt(disc) if(a.ne.0.0)then dnosemp = dnosemp/(2.0*a) dnosemp = dnosemp**0.33333 else dnosemp = 0.0 end if return end c ------------t end of nosedist t------------------------ c ----------t coord t---------------------------------------- subroutine coord(x,y,z,bx,by,bz,iswit) implicit none real*8 x,y,z,bx,by,bz integer iswit c c this subroutine transforms xm coord to GSM coord. c iswit=1 from model coord to GSM c iswit=-1 from GSM to model coord c if(iswit.eq.0) then x=-x y=-y endif c if(iswit.eq.1) then x=-x y=-y bx=-bx by=-by endif c if(iswit.eq.-1) then x=-x y=-y bx=-bx by=-by endif c return end c ----------t end of coord t--------------------------------- c -------------t rotat t---------------------------------- subroutine rotat(x,y,z,xr,yr,zr,rot) real*8 rot,x,y,z,xr,yr,zr,ct,st ct = cos(rot) st = sin(rot) xr = x yr = y*ct + z*st zr = z*ct - y*st return end c -------------t end of rotat t------------------------------- c --------------t atan3 t------------------------------------- real*8 function atan3(y,x) implicit double precision(a-h,o-z) atan3 = 0.00 if((x.ne.0.0).or.(y.ne.0.0)) atan3 = atan2(y,x) return end c ---------------------------------------------------------- c *********************************************************************** c subroutine norm (newer more efficient version ! 12/1/93 ) c this version only solves for a finite bn on the hemisphere c c subroutines include : norm, rotat, readnorm, readbessel, calcb, c dayn, legen, bcmpu, bcmpy, conv4, conv5 c *********************************************************************** c ------------t norm t------------------------------------------------ subroutine norm(stand,radius,iswit,x,y,z,bx,by,bz) implicit none integer ndim,nele,ndeg,ndegm,nstr,iswit real*8 stand,radius,xin,yin,zin,bx,by,bz,csph,snph,ph,ct,st, * bex,radius1,dt,pi,bxout,byout,bzout, x,y,z,rot parameter (pi=3.14159265358) common /store5/ csph,snph,ph common /store6/ ct,st common /store1/ radius1,bex,ndim,nele,ndeg,ndegm common /nparm/ dt,nstr save data ndim /99/,nele /2/ nstr=31 dt = 0.5/dfloat(nstr-1) nele = 1 c if iswit.eq.2 then calculate the b field if(iswit.eq.2)then call rotat(x,y,z,xin,yin,zin,rot) call calcb(xin,yin,zin,bxout,byout,bzout) call rotat(bxout,byout,bzout,bx,by,bz,-rot) return end if c if iswit.eq.1 then input the the data for COMMON /store1/ ndeg = 10 C write(6,*)' ndeg=',ndeg,' radius=',sngl(radius), C * ' nose distance =',sngl(stand) C write(6,*)' merging line rotation ',sngl(x) rot = x/180.0*pi radius1 = radius bex = radius - stand ndegm = min0(19,ndeg) if(iswit.eq.1)then C write(6,*) ' reading in the coefficients... ' C call readnorm C call readbessel return end if c otherwise, calculate the coefficients .............. c call begin c return write(6,'(a,$)')' Sorry, ' write(6,*)' this version cannot calculate the coefficients....' write(6,*)' future versions may do ' pause return end c -------------t end of norm t--------------------------- c ------------t readbessel t------------------------------------------ subroutine readbessel implicit double precision(a-h,o-z) common /bessel/ z0(20),ze(20,20),b0(20),be(20,20) open(unit=21,file='9bessel.dat',status='old') do 10 i=1,20 read(21,*) z0(i),b0(i) 10 continue do 11 i=1,20 do 12 j=1,20 read(21,*) ze(i,j),be(i,j) 12 continue 11 continue close(21) return end c ------------t end of readbessel t------------------------------- c -------------t readnorm t------------------------------ subroutine readnorm implicit none integer i,j,iin,jin real*8 g1in,g2in,axin,bxin,amin,bmin real*8 g1,g2,ax,bx,am,bm common /nosest/ g1(40,40),g2(40,40) common /store2/ ax(20,20),bx(20,20) common /zindep/ am(40),bm(40) open(unit=25,file='9norm.dat',status='old') do 11 i=1,40 do 11 j=1,i read(25,*)iin,jin,g1in,g2in if(iin.ne.i.or.jin.ne.j)then write(6,*)' !! Error reading in coeffs in "readnorm" ' stop end if g1(i,j) = g1in g2(i,j) = g2in 11 continue do 12 i=1,20 do 12 j=1,20 read(25,*)iin,jin,axin,bxin if(iin.ne.i.or.jin.ne.j)then write(6,*)' !! Error reading in coeffs in "readnorm" ' stop end if ax(i,j) = axin bx(i,j) = bxin 12 continue do 13 i=1,40 read(25,*)iin,amin,bmin if(iin.ne.i) then write(6,*)' !! Error reading in coeffs in "readnorm" !' stop end if am(i) = amin bm(i) = bmin 13 continue C write(6,*) ' ' close(25) return end c -------------t end of readnorm t------------------------------ c ------------t calcb t-------------------------------------------- subroutine calcb(xm,ym,zm,hxm,hym,hzm) implicit none integer ndim,nele,ndeg,ndegm real*8 xm,ym,zm,hxm,hym,hzm,csph,snph,phi, * pi,r,bex,ct,st,xv,yv,zv,d1,d2, * hxv,hyv,hzv,d3,rr,t,rrz,hr,hp,ht,hz parameter (pi=3.14159265358) common /store5/csph,snph,phi common /store1/ r,bex,ndim,nele,ndeg,ndegm common /store6/ ct, st call conv5(xm,ym,zm,xv,yv,zv,d1,d2,d3) call conv3(xv,yv,zv,rr,t) if(zv.ge.0.0) then rrz = dsqrt(xv*xv + yv*yv) call tailn(zv,rrz,hr,hz,hp) call bcmpy(phi,hr,hz,hp,hxv,hyv,hzv) else call dayn(t,rr,hr,ht,hp) call bcmpu(t,phi,hr,ht,hp,hxv,hyv,hzv) endif hxm = hzv hym = -hyv hzm = hxv return end c ------------t end of calcb t----------------------------- c ------------t dayn t--------------------------------------- c calculates the field in the sphere subroutine dayn(t,rr,hr,ht,hp) implicit none integer ndim,nele,ndeg,ndegm,m,n real*8 t,rr,hr,ht,hp,pnm,cpn,tpn,tpnm,spnm,gam1, * gam2,csph,snph,ph,c1,s1,pi,r,bex, * csphm,snphm,gama1,gama2,rrp,dhr,dhp,dht,rp,rp1, * dtot,fn,tiny,pn parameter (pi=3.14159265358,tiny=1.0d-15) common /leg1/pn(41),pnm(41,41),cpn(41),tpn(41),tpnm(41,41), * spnm(41,41) common /nosest/ gam1(40,40),gam2(40,40) common /store5/ csph,snph,ph common /store6/ c1, s1 common /store1/ r,bex,ndim,nele,ndeg,ndegm c pn = Legendre polynomials c cpn = derivative of pn by cos(theta) c tpn = derivative of pn by theta c pnm = modified Legendre polynomials c tpnm = derivative of pnm wrt cos(theta) c spnm = quotient pnm/sin(t) c hr - r component c ht = theta component c hp = phi component 1 call legen(t,ndeg,ndeg,2) if(rr.gt.r)go to 10 hr = (-gam1(1,1)*csph-gam2(1,1)*snph) * pnm(1,1) ht = (-gam1(1,1)*csph-gam2(1,1)*snph) * tpnm(1,1) hp = (gam1(1,1)*snph-gam2(1,1)*csph) * spnm(1,1) rp = 0.0d0 if (rr.gt.1.0d-3) rp = (rr/r) rrp = 1.0d0 do 3 n=2,ndeg rrp = rrp * rp fn = dfloat(n) do 4 m=1,n csphm = cos(dfloat(m)*ph) snphm = sin(dfloat(m)*ph) gama1= gam1(n,m) * rrp gama2= gam2(n,m) * rrp dhr=(-gama1*csphm-gama2*snphm)*pnm(n,m) dht=(-gama1*csphm-gama2*snphm)*tpnm(n,m)/fn dhp= (gama1*snphm-gama2*csphm)*spnm(n,m)/fn*m if(dabs(dhr).lt.tiny) dhr=0.0 if(dabs(dht).lt.tiny) dht=0.0 if(dabs(dhp).lt.tiny) dhp=0.0 dtot=dsqrt(dhr**2+(rr*dht)**2+(rr*s1*dhp)**2) if(dtot.lt.tiny) goto 5 hr = hr + dhr ht = ht + dht hp = hp + dhp 4 continue 3 continue 5 continue return 10 continue c ......... this part won't be called .................. rp1 =r/rr if (dabs(rp1).le.tiny) rp1=0 rp = rp1*rp1*rp1 gama1=gam1(1,1) gama2=gam2(1,1) hr = (-gama1*csph-gama2*snph)*pnm(1,1)*rp ht = (+gama1*csph+gama2*snph)*tpnm(1,1)*rp/2.0 hp = (-gama1*snph+gama2*csph)*spnm(1,1)*rp/2.0 do 20 n=2,ndeg rp =rp*rp1 fn=dfloat(n+1) do 20 m=1,n csphm=cos(ph*dfloat(m)) snphm=sin(ph*dfloat(m)) gama1=gam1(n,m)*rp gama2=gam2(n,m)*rp hr=hr+(-gama1*csphm-gama2*snphm)*pnm(n,m) ht=ht+(+gama1*csphm+gama2*snphm)*tpnm(n,m)/fn hp=hp+(-gama1*snphm+gama2*csphm)*spnm(n,m)/fn*m 20 continue 30 continue return end c ---------------t end of dayn t-------------------------- c ------------t legen t------------------------------------- c determines Legendre polynomial p-ne-me c for an angle t as well as the derivatives c and pnm/sin(t) subroutine legen(t,ne,me,imp) implicit double precision(a-h,o-z) parameter (pi=3.14159265358, tiny=1.0d-20) common /store6/ c1, s1 common /leg1/pn(41),pnm(41,41),cpn(41),tpn(41), * spnm(41,41),tpnm(41,41) common /store1/ r,bex,ndim,nele,ndeg,ndegm c pn = Legendre polynomials c cpn = derivatives of pn wrt cos(theta) c pnm = modified Legendre polynomials c tpnm = derivatives of Legendre polynomials wrt theta c spnm = pnm/sin(theta) c normalization based on Jahnke-Ende ****** nnnn = ne c1 = cos(t) s1 = sin(t) c2 = c1*c1 - s1*s1 s2 = 2.*s1*c1 pn(1) = c1 pn(2) = 0.25*(3.*c2+1.) cpn(1) = 1.0 cpn(2) = 3.*c1 tpn(1) = -s1 tpn(2) = -s1*cpn(2) pnm(1,1) = - tpn(1) sq3=sqrt(3.0) sq6=sqrt(6.0) sq2=sqrt(2.0) pnm(2,1) = + s1 * cpn(2)/sq3 pnm(2,2) = 3.0*s1*s1/(sq2*sq6) tpnm(2,2) = 6.0*s1*c1/(sq2*sq6) if(ne.lt.3) go to 20 if(ne.ge.3) nnnn= ne nnele = nnnn+1 do 1 n=3,nnele sqn = sqrt(n*(n+1.0)) zfe = 2.*n - 1.0 pn(n) = (1./n) * (zfe*c1*pn(n-1) - (n-1.0)*pn(n-2)) cpn(n) = zfe*pn(n-1) + cpn(n-2) tpn(n) = -s1 * cpn(n) pnm(n,1) = + s1*cpn(n)/sqn * sq2 if(abs(pnm(n,1)).le.tiny) pnm(n,1)=0.0 if(abs(cpn(n)).le.tiny) cpn(n)=0.0 if(abs(pn(n)).le.tiny) pn(n)=0.0 1 continue do 2 n = 1,nnnn sqn = sqrt(n*(n+1.0)) spnm(n,1)=cpn(n)/sqn * sq2 tpnm(n,1)=(sqn*pn(n) - c1*cpn(n)/sqn)*sq2 if(abs(spnm(n,1)).le.tiny) spnm(n,1)=0.0 if(abs(tpnm(n,1)).le.tiny) tpnm(n,1)=0.0 2 continue if(me.eq.1)return 20 nnele = nnnn+1 do 3 n=2,nnele mme=me+1 do 4 m=2,n sqnm=sqrt((n-m)*(n-m-1.0)) sqnm1=sqrt((n+m)*(n+m-1.0)) if(n-2.eq.0)go to 5 if(n-2.lt.m)go to 5 pnm(n,m) = pnm(n-2,m)*sqnm + (2.0*n-1.0)*s1*pnm(n-1,m-1) pnm(n,m)=pnm(n,m)/sqnm1 go to 40 5 pnm(n,m)= (2.0*n-1.0)*pnm(n-1,m-1)*s1/sqnm1 40 if(abs(pnm(n,m)).le.tiny) pnm(n,m)=0.0 4 continue 3 continue if(imp.eq.1)return do 7 n=2,ne do 8 m=2,n sqnm = sqrt((n+m)*(n-m+1.0)) sqnm1=sqrt((n-m)*(n+m+1.0)) if(n-1.lt.m)go to 9 tpnm(n,m)=0.5*(sqnm*pnm(n,m-1) - pnm(n,m+1)*sqnm1 ) go to 18 9 tpnm(n,m) = 0.5*(sqnm*pnm(n,m-1)) 18 sqnm = sqrt((n+m)*(n+m-1.0)) if(n-2.lt.m)go to 28 sqnm1 = sqrt((n-m)*(n-m-1.0)) spnm(n,m) = (sqnm*pnm(n-1,m-1)+pnm(n-1,m+1)*sqnm1)/(2.0*m) go to 80 28 spnm(n,m) = (sqnm*pnm(n-1,m-1))/(2.0*m) 80 if(abs(tpnm(n,m)).le.tiny) tpnm(n,m)=0.0 if(abs(spnm(n,m)).le.tiny) spnm(n,m)=0.0 8 continue 7 continue return end c -------------t end of legen t------------------- c ------------t bcmpu t--------------------------- subroutine bcmpu(t,phi,hrv,htv,hpv,hxv,hyv,hzv) implicit double precision(a-h,o-z) common /store5/ cp,sp,ph common /store6/ ct,st hxv = hrv*st*cp + htv*ct*cp - hpv*sp hyv = hrv*st*sp + htv*ct*sp + hpv*cp hzv = hrv*ct - htv*st return end c ----------------------------------------------------- c ----------------t bcmpy t------------------------------ subroutine bcmpy(phi,hrv,hzv,hpv,hxv,hyv,hzvv) implicit double precision(a-h,o-z) common /store5/ cp,sp,ph hxv = hrv*cp - hpv*sp hyv = hrv*sp + hpv*cp hzvv = hzv return end c ----------------------------------------------------- c --------------t conv3 t-------------------------------- subroutine conv3(xv,yv,zv,rr,t) implicit double precision(a-h,o-z) common /store5/ csph,snph,phi common /store6/ ct,st rr = sqrt(xv*xv + yv*yv + zv*zv) phi = atan3(yv,xv) t = atan3(sqrt(xv*xv+yv*yv),zv) ct = cos(t) st = sin(t) csph = cos(phi) snph = sin(phi) return end c ---------------t end of conv4 t---------------------- c ---------------t conv5 t----------------------------- subroutine conv5(xm,ym,zm,xv,yv,zv,xgse,ygse,zgse) implicit double precision(a-h,o-z) common /store1/ r,bex,ndim,nele,ndeg,ndegm xv = zm yv = -ym zv = xm - bex xgse = -xm ygse = -ym zgse = zm return end c ------------------------------------------------------------ c ---------------t tailn t---------------------------------------- c hhr = rho-omponent of the total field in the tail c hhp = phi-component of the total field in the tail c hhz = z-component of the total field in the tail subroutine tailn(z,rho,hhr,hhz,hhp) implicit none integer ndim,nele,ndeg,ndegm,i,mm real*8 z,rho,hhr,hhz,hhp,ax,bx,am,bm,bfjm2,bfjm1, * bfjm,z0,ze,b0,be,csph,snph,ph,pi,r,bex,ct,st, * d1,rrr,ex,ex1,ex2,b00,b11,b22,dhhr,dhhp, * dhhz,dto,arx,arz,cs,sn,bfj0,bfj1,bfj2,bfj,tiny parameter (pi=3.14159265358) common /store2/ ax(20,20),bx(20,20) common /zindep/ am(40),bm(40) common /bess/ bfjm2,bfjm1,bfjm common /bessel/ z0(20),ze(20,20),b0(20),be(20,20) common /store5/csph,snph,ph common /store1/ r,bex,ndim,nele,ndeg,ndegm common /store6/ ct,st data tiny /1.0d-10/, d1 /1.0d-6/ 50 continue hhr = 0.0 hhz = 0.0 hhp = 0.0 if(rho.gt.r)go to 120 c find the z independent solution c do 111 i=1,ndeg c cs = dcos(ph*i) c sn = dsin(ph*i) c if(i.eq.1)rrr=1. c if(i.ne.1)rrr=rrr*(rho/r) c if(dabs(am(i)).ge.d1)then c hhr = hhr+rrr*am(i)*cs c hhp = hhp-rrr*am(i)*sn c end if c if(dabs(bm(i)).ge.d1)then c hhr = hhr+rrr*bm(i)*sn c hhp = hhp+rrr*bm(i)*cs c end if c 111 continue c find the z dependent solution ...................... do 100 i = 1,ndegm cs = cos(ph*dfloat(i)) sn = sin(ph*dfloat(i)) if(i.eq.1)rrr=1. if(i.ne.1)rrr=rrr*(rho/r) do 101 mm = 1,ndegm arx = ze(i,mm)/r*rho arz = ze(i,mm)/r*z ex = 0.0d0 if(dabs(arz).ge.50.)go to 111 ex = exp(-arz) ex1 = ex * ax(i,mm)*ze(i,mm)/r ex2 = ex * bx(i,mm)*ze(i,mm)/r 111 continue b00 = 1.0d0 b11 = 0.0d0 b22 = 0.0d0 if(dabs(rho).lt.tiny) goto 112 if(i.gt.1) go to 113 b00 = bfj0(arx) b11 = bfj1(arx) b22 = bfj2(b11,b00,arx) go to 112 113 b22=bfj(i+1,arx) b11=bfjm1 b00=bfjm2 c ......... sign correction to 2nd term done 10/86 ............... 112 dhhr=-(b00-b22)/2.0*(ex1*cs + ex2*sn) dhhp= (b00+b22)/2.0*(ex1*sn - ex2*cs) dhhz= b11*(ex1*cs + ex2*sn) if(dabs(dhhr).lt.tiny) dhhr=0.0 if(dabs(dhhp*rho).lt.tiny) dhhp=0.0 if(dabs(dhhz).lt.tiny) dhhz=0.0 dto=dsqrt(dhhr**2+(dhhp*rho)**2+dhhz**2) if(dto.lt.tiny) go to 119 hhr = hhr + dhhr hhp = hhp + dhhp hhz = hhz + dhhz 101 continue 100 continue 119 continue return c find rho gt r solution .................... 120 do 121 i=1,40 cs = cos(ph*dfloat(i)) sn = sin(ph*dfloat(i)) if(i.eq.1) rrr=(r/rho)*(r/rho) if(i.ne.1) rrr=rrr*(r/rho) if(dabs(am(i)).lt.d1) goto 601 hhr = hhr+rrr*am(i)*cs hhp = hhp-rrr*am(i)*sn 601 if(dabs(bm(i)).lt.d1) goto 121 hhr = hhr+rrr*bm(i)*sn hhp = hhp+rrr*bm(i)*cs hhz = 0.0d0 121 continue return end c ---------------t end of tailn t------------------------------- c ----------t bfj t------------------------------------------------ c calculates Bessel functions c jn(x) by an iteration scheme real*8 function bfj(n,x) implicit double precision(a-h,o-z) common /bess/ bfj2,bfj11,bfjm d=1.0e-8 bfjm2=bfj0(x) bfj2=bfjm2 if(n.le.0)go to 1 bfjm1 = bfj1(x) bfj11=bfjm1 if(n.le.1)go to 2 if(abs(x).le.d)go to 3 if(n.eq.2)go to 77 c if((1.0*(n-2)/x)**(n-2).gt.1.0e04)go to 30 if((n-2.)*log((n-2.0)/x).gt.log(1.0e04))go to 30 77 do 10 k=2,n bfj=2.0*(k-1.0)/x*bfjm1 - bfjm2 bfj2=bfjm2 bfj11=bfjm1 bfjm2=bfjm1 bfjm1=bfj 10 continue if((1.0*n/x)**n.gt.1.0e04)bfj=0.0 if((1.0*(n-1)/x)**(n-1).gt.1.0e04)bfj11=0. if((n-1.)*log((n-1.0)/x).gt.log(1.0e04))bfj11=0. return 1 bfj=bfjm2 return 2 bfj=bfjm1 return 3 bfj = 0.0 return 30 bfj=0.0 bfj11=0.0 bfj2=0.0 return end c -----------t end of bfj t--------------------------- c -----------t bfj0 t---------------------------------- c ** j0(x) Bessel functions ***** c ** from Abramowitz + Stegun 9.4.1 and 9.4.3 ***** real*8 function bfj0(x) implicit double precision(a-h,o-z) data a1/-2.2499997/,a2/1.2656208/,a3/-0.3163866/,a4/0.0444479/, 1 a5/-0.0039444/,a6/0.00021/,b0/0.79788456/,b1/-0.00000077/, 2 b2/-0.0055274/,b3/-0.00009512/,b4/0.00137237/,b5/-0.00072805/, 3 b6/0.00014476/,c0/-0.78539816/,c1/-0.04166397/,c2/-0.00003954/, 4 c3/0.00262573/,c4/-0.00054125/,c5/-0.00029333/,c6/0.00013558/ x0=abs(x) t=x0/3.0 if (t.ge.1.0) goto 10 t=t*t h1=1.0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6))))) bfj0=h1 goto 20 10 t=1.0/t h2=b0+t*(b1+t*(b2+t*(b3+t*(b4+t*(b5+t*b6))))) h3=x0+c0+t*(c1+t*(c2+t*(c3+t*(c4+t*(c5+t*c6))))) bfj0=h2*cos(h3)/sqrt(x0) 20 continue return end c -----------t end of bfj0 t------------------------------ c -----------t bfj1 t------------------------------------- c ** j1(x) Bessel function ***** c ** from Abramowitz + Stegun 9.4.4 and 9.4.6 ***** real*8 function bfj1(x) implicit double precision(a-h,o-z) data a1/-.56249985/, 1 a2/.21093573/,a3/-.03954289/,a4/.00443319/,a5/-.00031761/, 2 a6/.00001109/,b0/.79788456/,b1/.00000156/,b2/.01659667/, 3 b3/.00017105/,b4/-.00249511/,b5/.00113653/,b6/-.00020033/, 4 c0/-2.35619449/,c1/.12499612/,c2/.0000565/,c3/-.00637879/, 5 c4/.00074348/,c5/.00079824/,c6/-.00029166/ x0=abs(x) rone = 1.0 x1=sign(rone,x) t=x0/3.0 if (t.ge.1.0) goto 10 t=t*t h1=0.5+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6))))) bfj1=x*h1 goto 20 10 t=1.0/t h2=b0+t*(b1+t*(b2+t*(b3+t*(b4+t*(b5+t*b6))))) h3=x0+c0+t*(c1+t*(c2+t*(c3+t*(c4+t*(c5+t*c6))))) bfj1=x1*h2*cos(h3)/sqrt(x0) 20 continue return end c -----------t end of bfj1 t------------------------------------- c -----------t bfj2 t--------------------------------------------- c ** j2(x) means the word b1 = j1(x) * * * * c ** b0 = j0(x) * * * * real*8 function bfj2(b1,b0,x) implicit double precision(a-h,o-z) bfj2 = 0.0 if(x.le.0.0) return bfj2 = (2./x) * b1 - b0 return end c -----------t end of bfj2 t--------------------------------------- c ************************************************************* c normp calculates the interconnection field in tail region c c subroutines include: normp, readinbn, calbn, rotat c ************************************************************* c ------------t normp t----------------------------------------- subroutine normp(stand,radius,xin,yin,zin, * bxout,byout,bzout,rat,iswit,ptl,ptl0) c implicit double precision(a-h, o-z) integer iswit,n real*8 xin,yin,zin, radius,stand,tl,x,y,z,xmp,ymp,zmp real*8 bx,by,bz,bxout,byout,bzout,bn, denom real*8 rho,rho2, cc, cc2, radius2, pi, rot, zf, tl2 real*8 tl0,tlv, rat, arg real*8 ptl,ptl0 parameter (pi=3.14159265358) parameter (n=5) save tl2 = radius bxout = 0.0 byout = 0.0 bzout = 0.0 bx = 0.0 by = 0.0 bz = 0.0 if(iswit.eq.1)then C write(6,*) 'Enter the information about the far tail xline : ' C write(6,*)' input tail length tl0,tlmp : (eg. 150, 300)' c read(5,*) tl,tl0 tl=ptl tl0=ptl0 C write(6,*)' tl0 =',sngl(tl),' tlmp =',sngl(tl0) C write(6,*)' tail xline rot =',sngl(xin) rot = xin/180.*pi c input bn field call readinbn return end if c calculation point call rotat(xin,yin,zin,x,y,z,rot) rho2 = radius**2-y**2 if(rho2.le.0.0) return rho = sqrt(rho2) cc = radius/10.0 cc2 = cc**2 radius2=radius**2 tlv = tl + (tl0-tl)*(y/radius)**2 rat = tlv/tl*(rho/radius)**n c ........... sphere region ........................... if(x.le.0.0) then c ........... inside magnetopuase retun ..................... if((x**2+y**2+z**2).le.radius**2) return if(x.le.-radius) return if(x**2+y**2.gt.radius**2) return c ........... only outside , non-zero ....................... xmp = x*sqrt((rho2+cc2)/(x**2+z**2+cc2)) ymp = y arg = rho2-xmp**2 if(arg.lt.0.0) return zmp = sqrt(rho2-xmp**2) c equation (3) th93 c bn dist varies as cos**(n+2) call calbn(xmp,ymp,zmp,bn) c ct = zmp/radius c if(ct.ne.0.0)bz=bn/ct c hyperbolic fit denom = rho2*(z**2+cc2)-x**2*cc2 if(denom.ne.0.0)then bx = bn*sqrt((rho2+cc2)/denom)* * x*z*rho/(x**2+z**2+cc2) bz = bn*sqrt((rho2+cc2)/denom)* * (z**2+cc2)/(x**2+z**2+cc2)*rho end if end if c ............ tail region ......................... if(x.gt.0.0) then xmp = x zf = 0.0 tlv = tl + (tl0-tl)*(y/radius)**2 if(rho.ne.0.0.and.x.lt.tlv) zf = rho*(1.0-x/(tlv)) if(abs(z).gt.zf)then if(abs(z).le.rho)then if(z.gt.0)then bx = - tlv/radius else bx = tlv/radius end if end if bz = (rho/radius) c **n dependence on bn see th93 equns 2a and 4 bx = bx*(rho/radius)**n bz = bz*(rho/radius)**n end if end if call rotat(bx,by,bz,bxout,byout,bzout,-rot) return end c --------------t end of normp t------------------------------ c --------------t calbn t---------------------------- c calculate the input bn subroutine calbn(xf,yf,zf,bnc) implicit double precision(a-h,o-z) common /block1/ bex,r,q,pi,psi common /bnor/ bn(31,121) pi = dacos(-1.0d0) c check to see if in the spherical regime, c if so, find the grid point; if not, only need the angle c from the symmetry plane... c i=1 is at the hemisphere/cylinder boundary c calculate the nearest grid point imax = 31 jmax = 121 c ............. if in spheric region ....................... if(xf.lt.0.)then c j point first if(yf.eq.0.0.and.zf.eq.0.0)then phi = 0.0 theta=pi/2.0 go to 11 end if phi = atan2(-yf,zf) if(phi.lt.0.)phi = phi+2.0*pi phi = (phi/(2.0*pi)*(jmax-1.0)) + 1. jphi = int(phi) c i point theta = asin((-xf)/r) 11 theta = (theta*2/pi*(imax-1.)) + 1. itheta = int(theta) c do 2-d interpolation using areas..... ithetp = itheta + 1 jphip = jphi + 1 djphip = float(jphip)-phi if(jphip.gt.jmax)jphip=1 if(itheta.eq.imax)then bnc = bn(imax,jmax) else area1 = (theta-float(itheta))*(phi-float(jphi)) area2 = (float(ithetp)-theta)*(phi-float(jphi)) area3 = (float(ithetp)-theta)*djphip area4 = (theta-float(itheta))*djphip if(area1.lt.0.0.or.area2.lt.0.0.or.area3.lt.0.0. * or.area4.lt.0.0)then write(6,*)' area lt 0' write(6,*)' area1 = ',area1 write(6,*)' area2 = ',area2 write(6,*)' area3 = ',area3 write(6,*)' area4 = ',area4 write(6,*)' x y z=',x,y,z end if bnc = area1*bn(ithetp,jphip) + * area2*bn(itheta,jphip) + * area3*bn(itheta,jphi) + * area4*bn(ithetp,jphi) end if c ............. if in cylinder region ....................... else c j point first again phi = atan2(-yf,zf) if(phi.lt.0.)phi = phi+2.0*pi phi = (phi/(2.0*pi)*(jmax-1.0)) + 1. jphi = int(phi) c do 1-d interpolation jphip = jphi + 1 djphip = float(jphip)-phi if(jphip.gt.jmax)jphip=1 bnc = bn(1,jphi)*(djphip) + * bn(1,jphip)*(phi-float(jphi)) end if return end c ------------t end of calbn t--------------------------------- c -------------t readinbn t------------------------------ c input the grids for dayside magnetopause subroutine readinbn implicit double precision(a-h,o-z) common /bnor/ bn(31,121) open(unit=22,file='9bnin.dat',status='old') do 120 i=1,121 read(22,22) (bn(j,i),j=1,31) 120 continue 22 format(3x,31f10.4) close(22) return end c ------------t end of readinbn t--------------------------------- c *********************************************************************** c subroutine tailst (last modified 2/86) c c subroutines include : tailst, rotat, geom, calct, g, h c *********************************************************************** c -----------t tailst t--------------------------------------- subroutine tailst(stand,radius,iswit,x,y,z,bx,by,bz) implicit double precision(a-h,o-z) parameter (pi=3.14159265358) common /store5/ csph,snph,ph common /store6/ ct, st common /store1/ r,bex,ndim,nele,ndeg,ndegm common /circ/ rd,c common /tailf/ bt,tw common /batch/ b,z11,z21,ib common /stretch/str,rlam,tlam,y0 r = radius bex = r - stand if(ib.eq.1) rot = tw/180.*3.14159 if(iswit.eq.2) then call rotat(x,y,z,xr,yr,zr,rot) call calct(xr,yr,zr,bxr,byr,bzr) call rotat(bxr,byr,bzr,bx,by,bz,-rot) return endif c if iswist=1 then c set up the geometry for the 'plasma sheet' c z1 is the thickness at y=0 and z2 at y=mp c rd is the radius of the fitted circle, c is the center bt = b z1 = z11 z2 = z21 if(ib.eq.1)then c write(6,*) ' tail field B =',sngl(b) c write(6,*) ' z11 =',sngl(z11),' z21 = ',sngl(z21) call geom(z11,z21) c write(6,*) ' rotation =',sngl(tw) return end if return end c -----------t end of tailst t----------------------------- c -----------t geom t--------------------------------------- subroutine geom(z1,z2) implicit double precision(a-h,o-z) parameter (pi=3.14159265358) common /store1/ r,bex,ndim,nele,ndeg,ndegm common /circ/ rd,zc rad=(r**2-z2**2) zc=(rad/(z2-z1)+(z1+z2))/2.0 rd=zc-z1 C write(6,*) ' Rd =',sngl(rd),' Zc =',sngl(zc) c parabolic fit c c=const, rd=slope c zc=z1 c rd=(z2-z1)/rad return end c --------------t end of geom t----------------------------- c --------------t calct t----------------------------------- subroutine calct(xm,ym,zm,bxm,bym,bzm) implicit double precision(a-h,o-z) parameter (pi=3.14159265358) common /store5/csph,snph,phi common /store1/r,bex,ndim,nele,ndeg,ndegm common /circ/ rd,cs common /store6/ ct, st common /stretch/str,rlam,tlam,y0 c .............. reset ..................................... bxm=0.0 bym=0.0 bzm=0.0 c determine whether inside magnetosphere .................. if((xm-bex).lt.0.0)then c ........... if dayside ................................... rrr=sqrt((xm-bex)**2+ym**2+zm**2) if(rrr.ge.r)return c th is a function of ym if(rd-abs(ym).lt.0.)then if(abs(ym).gt.r)then th=0.0 else th =sqrt(r**2-ym**2) end if else th=cs-sqrt(rd**2-ym**2) end if f=1.0 rr2=sqrt((xm-bex)**2+zm**2) if(rr2.lt.th) f = 1.0 +1.0*(rr2-th)/(th) if(rr2.ne.0.0)then bx = f*zm/rr2 by = 0.0 bz = -f*(xm-bex)/rr2 else bx = 0.0 by = 0.0 bz = 0.0 end if else c ........... nightside ...................................... rrz = sqrt(ym*ym + zm*zm) if(rrz.ge.r) return c uses functional dependence on ym gstr = g(str,rlam,y0,r,ym) c check if outside regions rrr=sqrt((gstr*(xm-bex))**2+(ym)**2+zm**2) if(rrr.ge.r)return c th is a function of ym if(rd-abs(ym).lt.0.)then if(abs(ym).gt.r)then th=0.0 else th =sqrt(r**2-(ym)**2) end if else th=cs-sqrt(rd**2-(ym)**2) end if f=1.0 uf = sqrt((gstr*(xm-bex))**2+zm**2) if(uf.lt.th) f = 1.0 +1.0*(uf-th)/(th) if(uf.ne.0.0)then bx = f*zm/uf by = 0.0 bz = -gstr*gstr*f*(xm-bex)/uf else bx = 0.0 by = 0.0 bz = 0.0 end if end if c y dependence on b factor = h(r,tlam,ym) bxm = bx * factor bym = by * factor bzm = bz * factor return end c -------------t end of calct t------------------------------- c ----------t g = y dependence of stretch t------------------- real*8 function g(str,rlam,y0,r,ym) implicit double precision(a-h,o-z) if(ym.lt.-y0) g = str*exp(-rlam*((ym+y0)/r)**2) if(ym.gt.y0) g = str*exp(-rlam*((ym-y0)/r)**2) if(abs(ym).le.y0) g = str return end c ------------------------------------------------------------ c -----------t h = y dependence of b t------------------------ real*8 function h(r,tlam,ym) implicit double precision(a-h,o-z) c parabola h = exp(-tlam*(ym/r)**2) return end c ---------------------------------------------------------- c ********************************************************** c mapping the electrical potential from the magnetopause c see Toffoletto and Hill 1989 paper. c c subroutines include : potin, rotat, calpot, calpotc c ********************************************************** c ---------------t potin t--------------------------------- subroutine potin(facin,angle) implicit double precision(a-h,o-z) parameter (pi=3.14159265358) common /pntl/ pot(31,121) common /calrot/ rot common /potrng/potmin,potmax open (unit=26,file='9potential.dat',status='old') do 120 i=1,121 read(26,222)(pot(j,i),j=1,31) 222 format(3x,31f10.4) 120 continue c normalize potmax = 0.0 potmin = 0.0 do 11 i=1,121 do 11 j=1,31 if(pot(j,i).ge.potmax)potmax=pot(j,i) if(pot(j,i).le.potmin)potmin=pot(j,i) 11 continue c write(6,*) ' ' c write(6,*) 'max pot =',sngl(potmax),' min pot =',sngl(potmin) c input potential rotation angle fac = facin*92.0 rot = angle C write(6,*)'Maximum polar cap potential drop in kV =',sngl(fac) C write(6,*)' potential rotation = ',sngl(rot) rot = rot/180.*pi c set to correct normalization do 12 i=1,121 do 12 j=1,31 pot(j,i) = - pot(j,i)/(potmax-potmin)*fac 12 continue close(26) c find min and max ranges potmax = pot(1,1) potmin = pot(1,1) do 200 i=1,31 do 200 j=1,121 if(potmax.le.pot(i,j))potmax=pot(i,j) if(potmin.ge.pot(i,j))potmin=pot(i,j) 200 continue C write(6,*) ' potmax=',sngl(potmax),' potmin=',sngl(potmin) C if(abs(potmax).gt.100.0.or.abs(potmin).gt.100.0) then C write(6,*) ' ' C write(6,*) ' AHA, YOU PUSH THE MODEL TOO FAR !' C write(6,*) ' the polar cap potential drop is > 100 KV !' C write(6,*) ' ' C pause C endif return end c ------------t end of potin t------------------------------ c ----------t calpot t------------------------------- subroutine calpot(xin,yin,zin,poten) implicit double precision(a-h,o-z) common /calrot/ rot common /block1/ bex,r,q,pi,psi common /pntl/ pot(31,121) c transform from GSM to model coordinates xin=-xin yin=-yin c check to see if you are in the spherical regime; c if so, find the grid point; if not, only need the angle c from the symmetry plane ... c rotate to specified angle call rotat(xin,yin,zin,xf,yf,zf,rot) c calculate the nearest grid point imax = 31 jmax = 121 if(xf-bex.lt.0.)then c j point first phi = atan2(-yf,zf) if(phi.lt.0.)phi = phi+2.0*pi if(phi.ge.2.0*pi)phi=phi-2.0*pi phi = phi/(2.0*pi)*dfloat(jmax-1) + 1. jphi = int(phi) c i point theta = asin((-xf+bex)/r) theta = theta*2.0/pi*dfloat(imax-1) + 1.0 itheta = int(theta) c do 2-d interpolation using areas..... ithetp = itheta + 1 jphip = jphi + 1 if(jphip.gt.jmax)jphip=1 if(itheta.eq.imax)then poten = pot(imax,jmax) else area1 = (theta-float(itheta))*(phi-float(jphi)) area2 = (float(ithetp)-theta)*(phi-float(jphi)) area3 = (float(ithetp)-theta)*(float(jphip)-phi) area4 = (theta-float(itheta))*(float(jphip)-phi) poten = area1*pot(ithetp,jphip) + * area2*pot(itheta,jphip) + * area3*pot(itheta,jphi) + * area4*pot(ithetp,jphi) end if else c j point first ...................................... phi = atan2(-yf,zf) if(phi.lt.0.0) phi = phi+2.0*pi if(phi.gt.2.0*pi) phi=phi-2.0*pi phi = phi/(2.0*pi)*dfloat(jmax-1) + 1.0 jphi = int(phi) c do 1-d interpolation ............................. jphip = jphi + 1 diff = float(jphip)-phi if(jphip.gt.jmax) jphip=1 poten = pot(1,jphi)*diff + * pot(1,jphip)*(phi-float(jphi)) end if c transform to GSM xin=-xin yin=-yin return end c ---------------t end of calpot t------------------------- c ---------------t calpotc t--------------------------------- subroutine calpotc(xin,yin,zin,poten) implicit double precision(a-h,o-z) common /potrng/potmin,potmax common /block1/ bex,r,q,pi,psi c transform from GSM to model coordinates xin=-xin yin=-yin poten = 0.0 rhoc = sqrt(yin**2+zin**2) c uses a simple approximation to find the potential as a c function of rho c for cos**3 dependence use - c poten = potmin*1.5*rhoc/r*(1.0-(rhoc/r)**2/3.) c poten = potmax*1.5*rhoc/r*(1.0-(rhoc/r)**2/3.) c for cos**5 dependence use - c poten = potmin*1.875*rhoc/r*(1.0-2.0*(rhoc/r)**2/3.+(rhoc/r)**4/5.) c poten = potmax*1.875*rhoc/r*(1.0-2.0*(rhoc/r)**2/3.+(rhoc/r)**4/5.) c for cos**7 dependence use - c poten=potmin*2.1875*rhoc/r c & *(1.0-(rhoc/r)**2+3.0/5.0*(rhoc/r)**4-1.0/7.0*(rhoc/r)**6) if(yin.le.0.0)then poten=potmin*2.1875*rhoc/r & *(1.0-(rhoc/r)**2+3.0/5.0*(rhoc/r)**4-(rhoc/r)**6/7.0) else poten=potmax*2.1875*rhoc/r & *(1.0-(rhoc/r)**2+3.0/5.0*(rhoc/r)**4-(rhoc/r)**6/7.0) end if c transform from GSM to model coordinates xin=-xin yin=-yin return end c ------------t end of calpotc t------------------------ c ********************************************************** c End of Toffoletto & Hill code c ********************************************************** c ********************************************************** c Calculate the magnetosheath B field c (taken from Spreiter & Stahara's code) c c Subroutines: init, bint, grid, quad c ********************************************************** c ---------s init s--------------------------------------------- c this subroutine initializes the array of the B-field with c the positions and the field vectors. subroutine init implicit double precision(a-h,o-z) common /imfb/ bimfx,bimfy,bimfz common /mpause/radius,stand,bex common /array/x(20,100),y(20,100),z(20,100) common /bfield/ bx(20,100),by(20,100),bz(20,100) common /bfieldx/ bxx(20,100),byx(20,100),bzx(20,100) common /parm/ dthet,dx,nblunt,nrad,nj,jm d = 3.14159265358/180. dthet = 4.0d0 nblunt = 24 nrad = 19 if(radius.le.0.0)then write(6,*)' !! Error in "init" : radius =< 0',radius stop end if open(unit=24,file='9msheathb.dat',status='old') bimfyz = sqrt(bimfy**2+bimfz**2) read(24,40)nrad,nj 40 format(2i3) if(bimfyz.ne.0.0)then do 1 j = 1,nj do 1 i = 1,nrad read(24,50)n,m,x(i,j),z(i,j),bx(i,j),bz(i,j),by(i,j) 50 format(2i3,5f10.4) x(i,j) = radius*x(i,j) z(i,j) = radius*z(i,j) bx(i,j) = -bimfyz*bx(i,j) bz(i,j) = -bimfyz*bz(i,j) by(i,j) = -bimfyz*by(i,j) 1 continue bx(1,1) = -bx(2,1) bz(1,1) = bz(2,1) do 5 j = 2,nj bmag1 = sqrt(bx(1,j)**2+bz(1,j)**2) bmag2 = sqrt(bx(2,j)**2+bz(2,j)**2) c renormalize i=1 case to i=2 case since I don't c believe the results for B at the magnetopause bx(1,j) = bx(1,j)/bmag1*bmag2 by(1,j) = by(2,j) bz(1,j) = bz(1,j)/bmag1*bmag2 5 continue end if c add bx contribution if(bimfx.ne.0.0)then do 11 j = 1,nj do 11 i = 1,nrad read(24,50) n,m,x0,z0,bxx(i,j),bzx(i,j),byx(i,j) bxx(i,j) = bimfx*bxx(i,j) byx(i,j) = 0.0 bzx(i,j) = bimfx*bzx(i,j) 11 continue end if close(24) return end c -------------s end of init s--------------------------------------- c -------------s bint s---------------------------------------------- c this subroutine determines the B-field c given the coordinates xp,zp c ***** note the sign of by has been changed*back!! ***** subroutine bint(xp,yp,zp,bbx,bby,bbz) implicit double precision(a-h,o-z) dimension z1(4),z2(4),f(4) common /imfb/ bimfx,bimfy,bimfz common /array/x(20,100),y(20,100),z(20,100) common /bfield/ bx(20,100),by(20,100),bz(20,100) common /bfieldx/ bxx(20,100),byx(20,100),bzx(20,100) common /parm/ dthet,dx,nblunt,nrad,nj,jm bbx = 0.0 bby = 0.0 bbz = 0.0 zpp = sqrt(yp**2+zp**2) call grid(xp,zpp,i,j) jm = j if(j.ge.89)return if(i.eq.20)go to 10 jinc = 1 jinc1 = 0 if(xp.ge.0.0)jinc=2 c if(j+1.eq.nblunt)jinc1 = 1 if(j+1.eq.nblunt)jinc= 2 if(j.le.2)jinc1 = 1 z1(3)=x(i,j+1-jinc) z1(4)=x(i+1,j) z1(2)=x(i,j+jinc1+jinc) z1(1)=x(i+1,j+1) z2(3)=z(i,j+1-jinc) z2(4)=z(i+1,j) z2(2)=z(i,j+jinc1+jinc) z2(1)=z(i+1,j+1) c do the extrapolation for the x-cpt f(3)=bx(i,j+1-jinc) f(4)=bx(i+1,j) f(2)=bx(i,j+jinc1+jinc) f(1)=bx(i+1,j+1) bbx = quad(z1,z2,f,xp,zpp) c do the extrapolation for the z-cpt f(3)=bz(i,j+1-jinc) f(4)=bz(i+1,j) f(2)=bz(i,j+jinc1+jinc) f(1)=bz(i+1,j+1) bbz = quad(z1,z2,f,xp,zpp) c do the extrapolation for the y-cpt f(3)=by(i,j+1-jinc) f(4)=by(i+1,j) f(2)=by(i,j+jinc1+jinc) f(1)=by(i+1,j+1) bby0= quad(z1,z2,f,xp,zpp) c do the extrapolation for the x-cpt from a Bx f(3)=bxx(i,j+1-jinc) f(4)=bxx(i+1,j) f(2)=bxx(i,j+jinc1+jinc) f(1)=bxx(i+1,j+1) bxxp= quad(z1,z2,f,xp,zpp) c do the extrapolation for the z-cpt from a Bx f(3)=bzx(i,j+1-jinc) f(4)=bzx(i+1,j) f(2)=bzx(i,j+jinc1+jinc) f(1)=bzx(i+1,j+1) bzxp= quad(z1,z2,f,xp,zpp) if(sqrt(yp**2+zp**2).ne.0.)then ang = atan2(yp,zp) else ang = 0.0d0 end if sina = sin(ang) cosa = cos(ang) bbx = bbx*cosa + bxxp bby =-(bby0 - bbz)*sina*cosa + bzxp*sina bbz = bbz*cosa**2 + bby0*sina**2 +bzxp*cosa return 10 bby = 0.0 bbz = -sqrt(bimfz**2+bimfy**2) bbx = -bimfx return end c ----------s end of bint s---------------------------------- c ------------s grid s---------------------------------------- c given an array and x,z coordinates this subroutine determines c between which grid points they lie. subroutine grid(xp,zp,i,j) implicit double precision(a-h,o-z) common /array/ x(20,100),y(20,100),z(20,100) common /parm/ dthet,dx,nblunt,nrad,nj,jm if(xp.gt.0.0)go to 50 j = nblunt if(abs(xp).le.1.0e-8) goto 11 j=1 axp = abs(xp) azp = abs(zp) ang = atan2(azp,axp)/3.14159*180. if(ang.le.dthet/2.0) goto 11 j = int((ang-dthet/2.0)/dthet+dthet/2.0) 11 do 15 i=1,nrad rr = sqrt(x(i,j)*x(i,j)+z(i,j)*z(i,j)) r = sqrt(xp*xp + zp*zp) if(r.le.rr) goto 20 15 continue i = nrad+2 20 i=i-1 if(i.eq.0)i=1 return 50 do 5 j = nblunt,nj if(abs(xp).le.abs(x(1,j))) goto 6 5 continue j = nj + 2 6 j = j-1 do 16 i=1,nrad if(abs(zp).le.abs(z(i,j)))go to 21 16 continue i = nrad+2 21 i=i-1 if(i.eq.0)i=1 return end c ----------s end of grid s----------------------------------- c ----------s quad s----------------------------------------- c this function performs a generalized quadrilateral interpolation, c given the coordinates of the vertices, the values of the function c at those points, and the coordinates of the point at which a c value of the function is desired, (xp,yp) c x and y are the coordinate arrays and f is the function array real*8 function quad(x,y,f,xp,yp) implicit double precision(a-h,o-z) dimension x(4),y(4),f(4) s1=4.0*xp-(x(1)+x(2)+x(3)+x(4)) s2=x(1)+x(2)-x(3)-x(4) s3=x(1)-x(2)+x(3)-x(4) s4=-x(1)+x(2)+x(3)-x(4) t1=4.0*yp-(y(1)+y(2)+y(3)+y(4)) t2=y(1)+y(2)-y(3)-y(4) t3=y(1)-y(2)+y(3)-y(4) t4=-y(1)+y(2)+y(3)-y(4) a=s2*t3-s3*t2 b=s1*t3+s2*t4-s3*t1-s4*t2 c=s1*t4-s4*t1 ds = b*b - 4.0*a*c d=sqrt(b*b-4.0*a*c) eta=-(b+d)/(a*2.0) if (abs(eta).gt.1.0) eta=-(b-d)/(2.0*a) xi=(s1+eta*s2)/(s4+eta*s3) quad=0.25*(f(1)*(1.0-eta)*(1.0-xi)+f(2)*(1.0-eta)*(1.0+xi) * +f(3)*(1.0+eta)*(1.0+xi)+f(4)*(1.0+eta)*(1.0-xi)) c for some reason this subroutine screws up sometimes giving absurd c values the problem seems to lie when eta is gt 1 and that b and d c are close as a safety default the following is inserted aquad=abs(quad) if(aquad.gt.abs(f(1)).and.aquad.gt.abs(f(2)).and. * aquad.gt.abs(f(3)).and.aquad.gt.abs(f(4)))quad=(f(1)+f(2) * +f(3)+f(4))/4.0 return end c ------------s end of quad s--------------------------------- c ********************************************************************* c the original name of this code, as programmed by G.-H. Voigt, c was "gehavo". this version, called "bvoigt", was modified by c R. V. Hilmer (8/87) to include a more complete set of comment c statements, including references to equation sources, etc. c the program's function is identical to the original version. c the main differences are cosmetic (i.e., subroutine name changes) c and are a consequence of translating the original german text c into english. c c for more information contact: G.-H. Voigt or R. V. Hilmer c Space physics and astronomy dept. c Rice University c P.O. Box 1892 c Houston, TX 77251 c c changes: 9/90 to subroutine cfdsph: increase filter size d-03 to d-02 c 11/90 to subroutine btata: Bessel filter size 170. to 160. c 02/94 converted to double precision and modified to work with c th93 c ********************************************************************* c c -------------v bvoigt v---------------------------------------------- subroutine bvoigt(stand,radius,tilt,cdipol,cimf,bimf,dir, 1 iswit, x,y,z, bx,by,bz) implicit double precision (a-h,o-z) c c purpose- computes three components of the vector magnetic field c at any point. points outside of the magnetopause will c only return nonzero fields if interconnection is allowed c with the imf. c c input- stand, the subsolar stand-off distance (in earth radii) c radius, radius of the magnetotail (in earth radii) c tilt, dipole tilt angle (in degrees) c tilt>0 indicates northern hemisphere summer c alfa, parameter indicating plasma content of tail, c = 1.0 means vacuum, = 0.0 means harris sheet c cdipol, fraction of dipole field allowed to penetrate c out through the magnetopause c cimf, fraction of the interplanetary field allowed to c penetrate in through the magnetopause c bimf, magnitude of the interplanetary magnetic field c (in gamma) c dir, direction of imf, 0 degrees is north, 90 degrees c is antisunward c iswit, indicates the procedure to follow on this call c bfield is called three times with iswit = 1,2 c then either 3 or 4 c *************altered so >2 gives gsm result 10/11/87 c = 1 general coefficients are found to represent c the magnetic field with legendre and bessel c series c = 2 constants specify the configuration desired c are calculated and stored, i.e. tilt c = 3 calculates field components, input and outpu c are in the sm (solar magnetospheric) system c = 4 same as 3 but gsm coordinates are used c x,y,z cartesian coordinate in either sm or gsm system c c output- bx,by,bz magnetic field vector in cartesian system c according to choice iswit=3 or 4. (in gamma) c c c variables- ndim number of grid points to be used for integration c to obtain series coefficients c nele number of iterations to be made between spherical c and cylindrical coefficient calculations. c bex distance between sphere center and earth c xd,yd,zd cartesian coordinates in dipole system, xd is c perpend. to dipole moment pointing antisunward c hxd,hyd,hzd magnetic field components in dipole system c xm,ym,zm earth centered cartesian system, xm points anti- c sunward, parallel to the tail axis (called the m c system) c hxm,hym,hzm magnetic field components on m system c snps,csps sin and cos of the dipole tilt angle c c references- c c ref1: voigt,g.-h., a mathematical magnetospheric c field model with independent physical parameters, planetary c and space science, vol. 29, 1-20,1981 c c ref2: abramowitz and stegun,eds., handbook of c mathematical functions, u.s. national bureau of standards, c dover, new york, (1965) c c ref3: jahnke-emde-losch, tables of higher c functions, 6th edition(revised), new york, (1960) c c ref4: voigt ph.d. thesis c c coordinate systems- c c gsm geocentric solar magnetospheric - earth centered c c xgsm : from earth to sun c ygsm : perpendicular to earth's dipole, pointing c towards dusk c zgsm : same sense as northern magnetic pole c dipole contained within x-z plane c c sm solar magnetic - earth centered c c xsm : toward sun but not directly unless tilt zero c ysm = ygsm c zsm : parallel to north magnetic pole c c m the m system (xm,ym,zm) - earth centered c c xm = -xgsm c ym = -ygsm c zm = zgsm c c v the v system (xv,yv,zv) - sphere centered c c xv = zgsm c yv = ygsm c zv = -xgsm - b where b = (tail radius)-(stand-off) c c d the d system (xd,yd,zd) - earth centered c c xd = -xsm c yd = -ysm = -ygsm c zd = zsm c c spherical system (r,theta,phi) - sphere centered c c r : distance from origin (xv,yv,zv)=(o,o,o) c theta : angle between zv and vector r c phi : angle between xv and component of vector r c in the xv-yv plane c c cylindrical system (rho,phi,z) - sphere centered c c rho : radial distance from zv axis c phi : angle between xv-zv plane and vector rho c (same as phi of spherical) c z = zv c common /block1/ bex,r,q,pi,psi common /block3/ h0, alfa, n common /block5/ csph,snph,qcs,qsn,csps,snps common /block6/ ct, st common /fimf/ himf, cd, ci, xi c r = radius himf = 0.0 xi = 0.0 cd = 0.0 cimf =0.0 c cd = cdipol c ci = cimf c himf = bimf c xi = dir c if (iswit.eq.1) then c get coefficients for spherical and bessel harmonics bex = radius - stand ndim = 101 nele = 4 call coeffs(ndim,nele) end if c if (iswit.eq.2) then c set constants, specify tilt which was preset in coeffs c note that sign of tilt angle is being switch ti = - tilt call const(ti) end if c c sm coordinate option cutoff for now everything in gsm ccc if (iswit.eq.3) then c get bfield at point using sm coordinate system c transform sm to dipole system, xd is perpendicular c to dipole moment pointing opposite xsm ccc xd = - x ccc yd = - y ccc zd = + z c c transform dipole to m system, xm antisolar ccc call trans4(xd,yd,zd,xm,ym,zm) c c get field components at xm,ym,zm ccc call btotal(xm,ym,zm,hxm,hym,hzm) c c transform field components back to dipole system ccc hxd = hxm*csps - hzm*snps ccc hyd = hym ccc hzd = hxm*snps + hzm*csps c c output in sm coordinates(gammas) ccc bx = - hxd ccc by = - hyd ccc bz = + hzd c ccc end if c c all calls with iswit > 2 give point result in gsm coordinates if (iswit.gt.2) then c get bfield at point using gsm coordinate system xm = - x ym = - y zm = + z c c get bfield components at xm,ym,zm call btotal(xm,ym,zm,hxm,hym,hzm) c c output in gsm coordinates(gammas) bx = - hxm by = - hym bz = + hzm c end if return end c -------------v end of bvoigt v--------------------------------------- c -------------v coeffs v---------------------------------------------- subroutine coeffs(ndim,nele) implicit double precision(a-h,o-z) c c purpose- calculate coefficients c input- ndim, number of grid points for numerical integration c nele, number of interations to obtain coefficients c c fill common block sites with bessel roots and values call vzero c do 10 ksa = 1, nele c get coefficients of legendre polynomials call sphrco(ndim,ksa) c get coefficients of bessel series call tailco(ndim) 10 continue c return end c -------------v end of coeffs v-------------------------------------- c -------------v const v---------------------------------------------- subroutine const(ti) implicit double precision(a-h,o-z) c c purpose- provide constants (to common blocks) c input- ti, dipole tilt angle in degrees c ti>0 indicates northern hemisphere winter c r, radius of tail (from block 1) in earth radii c output- q, ratio of dipole moment/(radius**3) in gamma c psi, dipole tilt angle in radians c pi, radian equivalent of 180 degrees c snps and csps, sin and cos of psi c variables- dm, dipole moment in gamma*(meters**3) c re, earth radius, meters c c common /block1/ bex,r,q,pi,psi common /block5/ csph,snph,qcs,qsn,csps,snps common /block6/ ct, st c pi = dacos(-1.00d0) c dm = 7.951d+24 re = 6.37104d+06 za = r*re q = dm/(za**3) c psi = (pi*ti)/180.0d0 csps = dcos(psi) snps = dsin(psi) c return end c -------------v end of const v---------------------------------------- c -------------v btotal v---------------------------------------------- subroutine btotal(xm,ym,zm,hxm,hym,hzm) implicit double precision(a-h,o-z) c c purpose- get field components using predetermined coefficients c input- xm,ym,zm m coordinates c output- hxm,hym,hzm field components in m system in gammas c common /block1/ bex,r,q,pi,psi common /block5/csph,snph,qcs,qsn,csps,snps common /block6/ ct, st c c transform m to spherical coordinates rr,t,phi c two step transform because zv is needed below in 'if' call trans5(xm,ym,zm,xv,yv,zv,d1,d2,d3) call trans3(xv,yv,zv,rr,t,phi) c if(zv.lt.0.0) then c get field from dayside of sphere zv<0 call daytot(rr,t,phi,hr,ht,hp) c transform to cartesian zv down tail call btran1(hr,ht,hp,hxv,hyv,hzv) else c get field from nightside of sphere zv>0 rrz = dsqrt(xv*xv + yv*yv) call tailto(rrz,phi,zv,hr,hp,hz) c transforn to cartesian zv down tail call btran2(hr,hz,hp,hxv,hyv,hzv) end if c c transform back to m system for output hxm = + hzv hym = - hyv hzm = + hxv c return end c -------------v end of btotal v--------------------------------------- c -------------v trans4 v---------------------------------------------- subroutine trans4(xd,yd,zd,xm,ym,zm) implicit double precision(a-h,o-z) c c purpose- transform dipole to m coordinates c input- xd,yd,zd dipole coordinates, zd north, xd tailward c output- xm,ym,zm m coordnates, xm is antisunward c variables- snps and csps, sin and cos of dipole tilt angle c common /block5/csph,snph,qcs,qsn,csps,snps c xm = zd*snps + xd*csps ym = yd zm = zd*csps - xd*snps c return end c -------------v end of trans4 v---------------------------------------------- c -------------v sphrco v---------------------------------------------- subroutine sphrco(ndim,ksa) implicit double precision(a-h,o-z) c c purpose- calculate coefficients for legendre polynomials c input- ndim, number of grid points for numerical integration c grid points are on sphere surface from theta 0 to c 180 and phi=0, theta=0 points down tail c ksa, iteration level c output- g1,g0, coefficients for tilt 0 and 90 degrees respect. c contain all except angular and legendre functions c see eq. 5.2.3 ref4 for general expressions c variables- q, ratio of dipole moment/radius**3 c r, radius of tail and sphere c b, = r - standoff distance c n, order of legendre polynomial c ti, dipole tilt angle c hs, grid size c y1,y0, integrand for 0 and 90 degree cases c z1,z0, result of integrals for 0 and 90 degree cases c common /block9/ x(1001),y1(1001),y0(1001),z1(1001),z0(1001), hs common /chablo/ g1(20), g0(20) common /block1/ b,r,q,pi,psi c c if(ksa.eq.1) then c coefficients are analytic, 5.11 and 5.12 ref1 ti = 0.0 call const(ti) g1(1) = - 2. * q g0(1) = - 2. * q ep = +1.0 rrp = 1.0 rp = b/r c do 1 n = 2, 20 ep = -ep rrp = rrp * rp g1(n) = - q * ep * float(n+1) * rrp g0(n) = - q * ep * float(n*(n+1)) * rrp 1 continue else c numerical integration required for ksa>1 c get radial component of dipole plus tail solution c on grid call brdip(ndim) do 3 n = 1, 20 c form integrands y1 and y0 call fun1(n,ndim) c integrate, results z1 and z0 call integr(hs,y1,z1,ndim) call integr(hs,y0,z0,ndim) zfe = float(2*n+1) fn = float(n) c form coefficients 5.2.3 ref4 g1(n) = 0.5 * zfe / (fn*(fn+1.)) * z1(ndim) g0(n) = 0.5 * zfe * z0(ndim) 3 continue end if c return end c -------------v end of sphrco v--------------------------------------- c -------------v daytot v---------------------------------------------- subroutine daytot(rr,t,phi,hr,ht,hp) implicit double precision(a-h,o-z) c c purpose- calculate total field of dayside half of sphere zv<0 c input- rr,t,phi spherical coordinates of interest c output- hr,ht,hp field components in spherical system c variables- cdipol, frac. of dipole penetrating out of magnetopause c cimf, frac. of imf penetrating into magnetopause c common /block1/ b, r, q, pi, psi common /fimf/ himf, cdipol, cimf, xi c cd = (1.- cdipol) cf = (1.- cimf) c if(rr.le.r) then c get field inside sphere c dipole field call bdip (rr,hr1,ht1,hp1) c chapman-ferraro due to dipole call cfdsph(rr,hr2,ht2,hp2) c tail solution, dayside influence call btasph (rr,hr3,ht3,hp3) c imf field c call bimf (himf,xi,d1,d2,d3,d4, hr5,ht5,hp5, d5,d6,d7) c c eq. 3.22 internal field ref1 c original version c hr = hr1 + cd * (hr2 + hr3) + cimf * hr5 c ht = ht1 + cd * (ht2 + ht3) + cimf * ht5 c hp = hp1 + cd * (hp2 + hp3) + cimf * hp5 hr = hr1 + cd * (hr2 + hr3) ht = ht1 + cd * (ht2 + ht3) hp = hp1 + cd * (hp2 + hp3) else c get field outside sphere c dipole field call bdip (rr,hr1,ht1,hp1) c chapman-ferraro due to imf c call cfisph(rr,hr2,ht2,hp2) c imf field c call bimf (himf,xi,d1,d2,d3,d4, hr5,ht5,hp5, d5,d6,d7) c c eq. 3.23 external field ref1 c original version c hr = hr5 + cf * hr2 + cdipol * hr1 c ht = ht5 + cf * ht2 + cdipol * ht1 c hp = hp5 + cf * hp2 + cdipol * hp1 hr = cdipol * hr1 ht = cdipol * ht1 hp = cdipol * hp1 end if return end c -------------v end of daytot v---------------------------------------- c -------------v tailto v---------------------------------------------- subroutine tailto(rr,phi,z,hr,hp,hz) implicit double precision(a-h,o-z) c c purpose- calculate total field of cylinder portion for zv>0 c input- rr,phi,z, cylindrical coordinates c output- hr,hp,hz, field components in cylindrical system c variables- cdipol, frac. of dipole penetrating out of magnetopause c cimf, frac. of imf penetrating in through magnetopause c snph,csph, sin and cos of phi c common /block1/ bex,r,q,pi,psi common /block5/ csph,snph,qcs,qsn,csps,snps common /block6/ ct, st common /fimf/ himf, cdipol, cimf, xi c cd = (1.- cdipol) cf = (1.- cimf) c if(rr.le.r) then c get field inside cylinder rrr = dsqrt(rr*rr + z*z) c dipole field call bdip(rrr,hr1,ht1,hp1) c transform spherical to cartesian, hzv down tail call btran1(hr1,ht1,hp1, hxv,hyv,hzv) c transform carteian to cylindrical hr1 = hxv * csph + hyv * snph hp1 = hyv * csph - hxv * snph hz1 = hzv c total tail field solution call btata(z,rr,hr2,hz2,hp2) c chapman-ferraro = total field - dipole field hr2 = hr2 - hr1 hp2 = hp2 - hp1 hz2 = hz2 - hz1 c get imf field c call bimf (himf,xi, d1,d2,d3,d4,d5,d6,d7, hr5,hz5,hp5) c sum of internal fields, 3.22 ref1 c orignal version c hr = hr1 + cd * hr2 + cimf * hr5 c hz = hz1 + cd * hz2 + cimf * hz5 c hp = hp1 + cd * hp2 + cimf * hp5 hr = hr1 + cd * hr2 hz = hz1 + cd * hz2 hp = hp1 + cd * hp2 else c get field outside of cylinder rrr = dsqrt(rr*rr + z*z) c dipole field call bdip(rrr,hr1,ht1,hp1) c transform spherical to cartesian, hzv down tail call btran1(hr1,ht1,hp1, hxv,hyv,hzv) c transform cartesian to cylindrical hr1 = hxv * csph + hyv * snph hp1 = hyv * csph - hxv * snph hz1 = hzv c chapman-ferraro due to imf c call cfita(rr,hr2,hz2,hp2) c imf field c call bimf (himf,xi, d1,d2,d3,d4,d5,d6,d7, hr5,hz5,hp5) c sum of external fields using 3.23 ref1 c hr = hr5 + cf * hr2 + cdipol * hr1 c hz = hz5 + cf * hz2 + cdipol * hz1 c hp = hp5 + cf * hp2 + cdipol * hp1 hr = cdipol * hr1 hz = cdipol * hz1 hp = cdipol * hp1 end if return end c -------------v end of tailto v--------------------------------------- c -------------v btran1 v---------------------------------------------- subroutine btran1(hrv,htv,hpv,hxv,hyv,hzv) implicit double precision(a-h,o-z) c c purpose- changes field components from spherical to cartesian c hzv points down tail axis, =0 to sphere center c input- hrv,htv,hpv, spherical components c output- hxv,hyv,hzv, cartesian components c variables- st,ct, sin and cos of theta c cp,sp, cos and sin of phi c common /block5/ cp,sp,qcs,qsn,csps,snps common /block6/ ct, st c hxv = hrv*st*cp + htv*ct*cp - hpv*sp hyv = hrv*st*sp + htv*ct*sp + hpv*cp hzv = hrv*ct - htv*st c return end c -------------v end of btran1 v--------------------------------------- c -------------v tailco v---------------------------------------------- subroutine tailco(ndim) implicit double precision(a-h,o-z) c c purpose- numerical calculation of fourier-bessel coefficients c input- ndim, number of grid points along cylinder cap at c theta=90 degrees for integration c output- ax,bx fourier-bessel coefficients with modifications c to be useful in final expression 6.4 ref1 c variables- be,bo bessel functions j1 and j0 with arguments that c are roots of derivatives of the bessel functions c see 5.21 and 5.22 ref1 c z1,z0 results of integration for 0 and 90 degrees cases c ze,zo zeros of derivatives of j1 and j0 c q, dipole moment/(radius**3) c common /block9/ x(1001),y1(1001),y0(1001),z1(1001),z0(1001), hs common /block1/ b ,r, q, pi, psi common /block2/ ze(20), ax(20), zo(20), bx(20) common /block7/ be(20), bo(20) c c generate potential on cylinder head= chap-ferr + dipole call tapot(ndim) c generate coefficients for n orders do 1 n = 1, 20 c form integrands y1 and y0 call fun2(n,ndim) c do integrations for both tilt cases,result z1, z0 call integr(hs,y1,z1,ndim) call integr(hs,y0,z0,ndim) c prepare some bessel factors bx0 = be(n) * be(n) by0 = bo(n) * bo(n) zex = ze(n) * ze(n) c coefficients a and b of 5.24 and 5.25 ref1 ax(n) = z1(ndim) * zex / (bx0 * (zex-1.)) bx(n) = z0(ndim) / (by0) c modify coefficients to suit expression 6.4 ref1 ax(n) = ax(n) * ze(n) * q bx(n) = bx(n) * zo(n) * 2. * q 1 continue c return end c -------------v end of tailco v--------------------------------------- c -------------v trans3 v---------------------------------------------- subroutine trans3(xv,yv,zv,rr,t,phi) implicit double precision(a-h,o-z) c c purpose- change cartesian to spherical and form trig. functions c for common blocks c input- xv,yv,zv cartesian with zv down tail c output- rr,t,phi spherical coordinates c variables- q, ratio dipole moment/(radius**3) c common /block1/ bex,r,q,pi,psi common /block5/ csph,snph,qcs,qsn,csps,snps common /block6/ ct, st c rr = dsqrt(xv*xv + yv*yv + zv*zv) phi = vatan3(yv,xv) t = vatan3(dsqrt(xv*xv+yv*yv),zv) c ct = dcos(t) st = dsin(t) csph = dcos(phi) snph = dsin(phi) qcs = q*csph qsn = q*snph c return end c -------------v end of trans3 v--------------------------------------- c -------------v trans5 v---------------------------------------------- subroutine trans5(xm,ym,zm, xv,yv,zv, xgse,ygse,zgse) implicit double precision(a-h,o-z) c c purpose- transform m to spherical and gse system c input- xm,ym,zm the m system, earth centered, xm down tail c output- xv,yv,zv cartesian, sphere centered, zv down tail c xgse,ygse,zgse, gse system c variables- bex, = radius - stand-off distance c common /block1/ bex,r,q,pi,psi xv = zm yv = -ym zv = xm - bex c xgse = - xm ygse = -ym zgse = + zm c return end c -------------v end of trans5 v-------------------------------------- c -------------v vzero v---------------------------------------------- subroutine vzero implicit double precision(a-h,o-z) c c purpose- fill common blocks with zeros of bessels and function c values at zeros c input- source is table 9.5 ref2 c ze,zeros of j1'(xi) c zo, zeros of jo'(yi) c be, bessel value j1(xi) c bo, besel value jo(yi) c common /block2/ ze(20),ax(20),zo(20),bx(20) common /block7/ be(20),bo(20) c c c zeros j1'(xi) = 0.0 c ze(1)= 1.84118 ze(2)= 5.33144 ze(3)= 8.53632 ze(4)= 11.70600 ze(5)= 14.86359 ze(6)= 18.01553 ze(7)= 21.16437 ze(8)= 24.31133 ze(9)= 27.45705 ze(10)= 30.60192 ze(11)= 33.74618 ze(12)= 36.88999 ze(13)= 40.03344 ze(14)= 43.17663 ze(15)= 46.31960 ze(16)= 49.46239 ze(17)= 52.60504 ze(18)= 55.74757 ze(19)= 58.89000 ze(20)= 62.03235 c c c zeros j0'(yi) = 0.0 c zo(1) = 3.8317059702 zo(2) = 7.0155866698 zo(3) = 10.1734681351 zo(4) = 13.3236919363 zo(5) = 16.4706300509 zo(6) = 19.6158585105 zo(7) = 22.7600843806 zo(8) = 25.9036720876 zo(9) = 29.0468285349 zo(10) = 32.1896799110 zo(11) = 35.3323075501 zo(12) = 38.4747662348 zo(13) = 41.6170942128 zo(14) = 44.7593189977 zo(15) = 47.9014608872 zo(16) = 51.0435351836 zo(17) = 54.1855536411 zo(18) = 57.3275254379 zo(19) = 60.4694578453 zo(20) = 63.61141 c c c be(i) = j1(xi) bessels evaluated at xi's of above c be(1)= +0.58187 be(2)= -0.34613 be(3)= +0.27330 be(4)= -0.23330 be(5)= +0.20701 be(6)= -0.18802 be(7)= +0.17346 be(8)= -0.16184 be(9)= +0.15228 be(10)= -0.14424 be(11)= +0.13736 be(12)= -0.13137 be(13)= +0.12611 be(14)= -0.12143 be(15)= +0.11724 be(16)= -0.11345 be(17)= +0.11001 be(18)= -0.10687 be(19)= +0.10397 be(20)= -0.10131 c c c bo(i) = j0(yi) bessels evaluated at yi's of above c bo(1) = -0.4027593957 bo(2) = +0.3001157525 bo(3) = -0.2497048871 bo(4) = +0.2183594072 bo(5) = -0.1964653715 bo(6) = +0.1800633753 bo(7) = -0.1671846005 bo(8) = +0.1567249836 bo(9) = -0.1480111100 bo(10) = +0.1406057982 bo(11) = -0.1342112403 bo(12) = +0.1286166221 bo(13) = -0.1236679608 bo(14) = +0.1192498120 bo(15) = -0.1152736941 bo(16) = +0.1116704969 bo(17) = -0.1083853489 bo(18) = +0.1053740554 bo(19) = -0.1026005671 bo(20) = +0.1000 c c return end c -------------v end of vzero v------------------------------------------ c -------------v bimf v-------------------------------------------------- subroutine bimf(himf,xi, hxm,hzm, hxv,hzv, hr,ht,hp, hrho,hz,hphi) implicit double precision(a-h,o-z) c c purpose- determine imf contribution in cartesian m and v,spher, c and cylindrical coordinates c no y component included c input- himf, magnitude of imf (gamma) c xi, direction of imf, 0 degrees is northward, 90 antisun c output- hxm,hzm field values in m system, hxm down tail c hxv,hzv field values in v system, hzv down tail c hr,ht,hp " " " spherical c hrho,hz,hphi " " cylindrical c variables- st,ct sin and cos of theta c snph,csph sin and cos of phi c common /block5/csph,snph,qcs,qsn,csps,snps common /block6/ ct, st c degrees to radians axi = xi*3.1415926535/180.0 c m system hxm = himf * dsin(axi) hzm = himf * dcos(axi) c v system hxv = hzm hzv = hxm c spherical system hr = + hxv*st*csph + hzv*ct ht = + hxv*ct*csph - hzv*st hp = - hxv*snph c cylindrical system hrho = + hxv * csph hphi = + hp hz = + hzv c return end c ----------v end of bimf v--------------------------------------------- c ----------v cfdsph v-------------------------------------------------- subroutine cfdsph(rr,hr,ht,hp) implicit double precision(a-h,o-z) c c purpose- calculate field of chapman-ferraro currents for c shielding dipole only in sphere. summation of field due c to potential 5.11 ref1 c note: recurrence relations Jahnke-Emde normalized, ref3 c input- rr, radial location c output- hr,ht,hp spherical field components, 5.2.13 ref4 c variables- gam1,gamo legendre polynomial coefficients from sphrco c pn, legendre polynomial m=0 c cpn, derivative of pn w.r.t. argument dcos(theta) c tpn, derivative of pn w.r.t. theta c pnm, legendre polynomial m=1 c tpnm, derivative of pnm w.r.t. theta c spnm, quotient pn/dsin(theta) c s1,c1 sin and cos of theta c dimension pn(20), cpn(20), tpn(20) dimension pnm(20),tpnm(20),spnm(20) common /chablo/ gam1(20), gam0(20) c common /block1/ b,r,q,pi,psi common /block5/csph,snph,qcs,qsn,csps,snps common /block6/ c1, s1 c d = 1.0d-02 c c2 = c1*c1 - s1*s1 s2 = 2.*s1*c1 c pn(1) = c1 pn(2) = 0.25*(3.*c2+1.) cpn(1) = 1.0 cpn(2) = 3.*c1 c do 1 n=3,20 fn = float(n) zfe = 2.*fn - 1.0 pn(n) = (1./fn) * (zfe*c1*pn(n-1) - (fn-1.)*pn(n-2)) cpn(n) = zfe*pn(n-1) + cpn(n-2) 1 continue c do 2 n=1,20 fn = float(n) tpn(n) = - s1*cpn(n) pnm(n) = - tpn(n) spnm(n) = + cpn(n) tpnm(n) = fn*(fn+1.)*pn(n) - c1*cpn(n) 2 continue c hr1 = gam1(1) * pnm(1) hr2 = gam0(1) * pn(1) ht1 = gam1(1) * tpnm(1) ht2 = gam0(1) * tpn(1) hp = gam1(1) * spnm(1) c rp = 0.0 if(rr.gt.d) rp = (rr/r) rrp = 1.0 c do 3 n=2,20 fn = float(n) c rrp = rrp * rp c gama = gam1(n) * rrp gamb = gam0(n) * rrp c c 1 and 2 indicate terms meant for 0 and 90 degree cases hr1 = hr1 + gama * pnm(n) hr2 = hr2 + gamb * pn(n) c ht1 = ht1 + gama * tpnm(n)/ fn ht2 = ht2 + gamb * tpn(n) / fn c hp = hp + gama * spnm(n)/ fn c 3 continue c c each part relating to tilt cases is combined with c functions of 5.9 ref1 hr = - (hr1 * csph * csps + hr2 * snps) ht = - (ht1 * csph * csps + ht2 * snps) hp = + hp * snph * csps c c return end c ----------v end of cfdsph v--------------------------------------------- c ----------v cfisph v---------------------------------------------------- subroutine cfisph(rr,hr,ht,hp) implicit double precision(a-h,o-z) c c purpose- calculate chapman-ferraro field outside shere due to imf c using 5.15 ref1 c input- rr, radial coordinate c output- hr,ht,hp chap-ferr field components in spherical system c variables- hr3,ht3,hp3 imf field components from subrout. bimf c c common /block1/ bex, r, q, pi, psi common /fimf/ himf, cdipol, cimf, xi c c get imf field components call bimf(himf,xi, d1,d2,d3,d4, hr3,ht3,hp3, d5,d6,d7) c transform components to spherical system fr = r/rr fr = fr*fr*fr c hr = - hr3 * fr ht = + 0.5 * ht3 * fr hp = + 0.5 * hp3 * fr c return end c --------------v end of cfisph v---------------------------------------- c --------------v bdip v------------------------------------------------- subroutine bdip(rr,hr,ht,hp) implicit double precision(a-h,o-z) c c purpose- calculate explicit components of dipole using spherical c equivalent of 5.6 ref1 for potential c input- rr, radial coordinate location c output- hr,ht,hp dipole field in spherical system (gamma) c variables- c1,s1 cos and sin of theta c r, radius of sphere and tail c b, = r - standoff distance c snps,csps sin and cos of psi, the dipole tilt angle c common /block1/ b,r,q,pi,psi common /block5/csph,snph,qcs,qsn,csps,snps common /block6/ c1, s1 c r3 = r**3 c a = qcs*r3 aa = qsn*r3 c v = ((rr*rr) + (b*b) + (2.*b*rr*c1)) c av = dabs(v) sv = dsign(1.0d0,v) c w3 = dsqrt(av**3) w4 = dsqrt(av**5) c w1 = w3 * sv w2 = w4 * sv c f1 = 2.*rr*b*s1 g1 = 2.*rr + 2.*b*c1 e1 = 1.5*rr*g1/w2 e2 = 1.5*f1/w2 c hr1 = a*csps*s1*(1./w1 - e1) hr2 = q*r3*c1*snps*(1./w1 - e1) hr3 = -q*r3*b*snps*3.*(rr+ b*c1)/w2 c ht1 = a*csps*(c1/w1 + s1*e2) ht2 = q*r3*snps*(c1*e2 - s1/w1) ht3 = q*r3*b*snps*(3.*b*s1/w2) c hp = -aa*csps/w1 c hr = hr1 + hr2 + hr3 ht = ht1 + ht2 + ht3 hp = hp c return end c --------------v end of bdip v-------------------------------------- c --------------v fun1 v--------------------------------------------- subroutine fun1(n,ndim) implicit double precision(a-h,o-z) c c purpose- form integrand for coefficient calculations at each of c ndim steps for a certain order of legendre polynomial c input- n, order of legendre polynomial c ndim, number of grid points used for integration c output- y1,yo value of integrand for cases 0 and 90 degrees resp c variables- pnm, legendre polynomial m=1 c pn, legendre polynomial m=0 c st, dsin(t) last value comes from subrout. legen c common /vleg1/ pn(20), pnm(20), cpn(20) common /block8/ w1(1001),w0(1001) common /block9/ x(1001),y1(1001),y0(1001),z1(1001),z0(1001), hs common /block6/ ct, st c c do 1 i = 1, ndim c x(i) is theta, generates legendre polynomial terms call vlegen(x(i),n) c d = 1.0d-20 c if(abs(w1(i)).lt.d) w1(i) = 0.0 if(abs(w0(i)).lt.d) w0(i) = 0.0 if(abs(pnm(n)).lt.d) pnm(n) = 0.0 if(abs(pn(n)).lt.d) pn(n) = 0.0 if(abs(st).lt.d) st = 0.0 c if(abs(pnm(n)*w1(i)).lt.d) pnm(n) = 0.0 if(abs(pn(n)*w0(i)).lt.d) pn(n) = 0.0 c c ndim integrands for each order to get coefficients c 5.2.3 ref4 y1(i) = w1(i) * pnm(n) * st y0(i) = w0(i) * pn(n) * st c 1 continue c c return end c -----------v end of fun1 v--------------------------------------- c -----------v brdip v--------------------------------------------- subroutine brdip(ndim) implicit double precision(a-h,o-z) c c purpose- generate radial components of dipole field on grid for c tilt 0 and 90 degrees c input- ndim, number of grid points at radius=r along curve c theta=0 to 180 degrees to be used for integration c output- w1,wo radial components of dipole, into block 8 c variables- kdim, number of the grid point separating day and night- c side of sphere c r, radius of shpere and tail c q, ratio dipole moment/(r**3) c snph,csph sin and cos of phi c qsn,qcs =q*snph and =q*csph respectively c x(i), theta values of grid points c hs, grid point spacing c common /block1/ b,r,q,pi,psi common /block3/ h0,alfa,kswit common /block5/ csph,snph,qcs,qsn,csps,snps common /block6/ ct,st common /block8/ w1(1001),w0(1001) common /block9/ x(1001),y1(1001),y0(1001),z1(1001),z0(1001), hs c kdim = (ndim-1)/2 c c phi = 0 degrees for this subrout. csph = 1.0 snph = 0.0 qcs = q qsn = 0.0 c c c generate grid points x(1) = 0.0 hs = pi/float(ndim-1) do 1 i = 2, ndim x(i) = float(i-1) * hs 1 continue c c get components for tilt=0 degrees ti = 0.0 call const(ti) c do 4 i = 1, ndim c ct = dcos(x(i)) st = dsin(x(i)) c rrdip = r c get radial component of dipole at sphere radius call bdip(rrdip,hr,ht,hp) w1(i) = hr c if (i.le.kdim) then c location is nightside, change to cylindrical rho= r * st z = r * ct c get cylinder solution call vacta(z,rho,hr,hz,hp) hr = hr * st + hz * ct c cylinder solution has dipole + chap-ferr 5.16 ref1 w1(i) = w1(i) - hr else c location is dayside of sphere rrshe = r c dayside field resulting from tail solution call btasph(rrshe,hr,ht,hp) c dipole plus dayside tail solution gives dayside c field, 5.1 ref1 w1(i) = w1(i) + hr end if 4 continue c c c generate components for tilt=90 degrees c (same procedure as above for tilt=0 degrees) ti = 90.0 call const(ti) c c do 7 i = 1, ndim c ct = dcos(x(i)) st = dsin(x(i)) c rrdip = r call bdip(rrdip,hr,ht,hp) w0(i) = hr c if (i.le.kdim) then rho= r * st z = r * ct call vacta(z,rho,hr,hz,hp) hr = hr * st + hz * ct w0(i) = w0(i) - hr else rrshe = r call btasph(rrshe,hr,ht,hp) w0(i) = w0(i) + hr end if 7 continue c c return end c ------------v end of brdip v--------------------------------------- c ------------v integr v--------------------------------------------- subroutine integr(h,y,z,ndim) implicit double precision(a-h,o-z) c c purpose- general integration routine c input- h, grid size c y, integrand, from fun1 or fun2 subrout. c ndim, number of grid points used c output- z integral result of function y from zero to ndim*h c (used to fill z1 and zo of block 9) c dimension y(ndim), z(ndim) c ht=.3333333*h if(ndim-5)7,8,1 c c ndim is greater than 5. preparations of integration loop 1 sum1=y(2)+y(2) sum1=sum1+sum1 sum1=ht*(y(1)+sum1+y(3)) aux1=y(4)+y(4) aux1=aux1+aux1 aux1=sum1+ht*(y(3)+aux1+y(5)) aux2=ht*(y(1)+3.875*(y(2)+y(5))+2.625*(y(3)+y(4))+y(6)) sum2=y(5)+y(5) sum2=sum2+sum2 sum2=aux2-ht*(y(4)+sum2+y(6)) z(1)=0. aux=y(3)+y(3) aux=aux+aux z(2)=sum2-ht*(y(2)+aux+y(4)) z(3)=sum1 z(4)=sum2 if(ndim-6)5,5,2 c c integration loop 2 do 4 i=7,ndim,2 sum1=aux1 sum2=aux2 aux1=y(i-1)+y(i-1) aux1=aux1+aux1 aux1=sum1+ht*(y(i-2)+aux1+y(i)) z(i-2)=sum1 if(i-ndim)3,6,6 3 aux2=y(i)+y(i) aux2=aux2+aux2 aux2=sum2+ht*(y(i-1)+aux2+y(i+1)) 4 z(i-1)=sum2 5 z(ndim-1)=aux1 z(ndim)=aux2 return 6 z(ndim-1)=sum2 z(ndim)=aux1 return c end of integration loop c 7 if(ndim-3)12,11,8 c c ndim is equal to 4 or 5 8 sum2=1.125*ht*(y(1)+y(2)+y(2)+y(2)+y(3)+y(3)+y(3)+y(4)) sum1=y(2)+y(2) sum1=sum1+sum1 sum1=ht*(y(1)+sum1+y(3)) z(1)=0. aux1=y(3)+y(3) aux1=aux1+aux1 z(2)=sum2-ht*(y(2)+aux1+y(4)) if(ndim-5)10,9,9 9 aux1=y(4)+y(4) aux1=aux1+aux1 z(5)=sum1+ht*(y(3)+aux1+y(5)) 10 z(3)=sum1 z(4)=sum2 return c c ndim is equal to 3 11 sum1=ht*(1.25*y(1)+y(2)+y(2)-.25*y(3)) sum2=y(2)+y(2) sum2=sum2+sum2 z(3)=ht*(y(1)+sum2+y(3)) z(1)=0. z(2)=sum1 12 return end c --------------v end of integr v---------------------------------------- c -------------v btran2 v------------------------------------------------ subroutine btran2(hrv,hzv,hpv,hxv,hyv,hzvv) implicit double precision(a-h,o-z) c c purpose- change field components from cylindrical to cartesian c input- hrv,hzv,hpv cylindrical field components, hzv down tail c output- hxv,hyv,hzvv cartesian components, hzvv down tail c variables- sp,cp sin and cos of phi c common /block5/ cp,sp,qcs,qsn,csps,snps c hxv = hrv*cp - hpv*sp hyv = hrv*sp + hpv*cp hzvv = hzv c return end c ------------v end of btran2 v------------------------------------------ c ------------v btasph v------------------------------------------------- subroutine btasph(rr,hr,ht,hp) implicit double precision(a-h,o-z) c c purpose- calculate field components of tailfield within sphere c for zv<0. (image of vacta) components derived from 5.26 c ref1 with sign change in exponent due to zv position c first called from sphrco on iteration 2. 15 terms used. c input- rr, radial distance in spherical coordinates c output- hr,ht,hp spherical field components c variables- rho,z cylindrical equivalents of input c st,ct sin and cos of theta c ze, zeros of j1'( ) bessel 5.23 ref1 c zo, zeros of jo'( ) bessel " c ax,bx fourier-bessel coefficients with additions c alpha, plasma content of tail 1 means vacuum c 0 means harris sheet c alpha=1 used to get coefficients c alpha not = 1 only possible for btotal application c c common /block1/ b,r,q,pi,psi common /block2/ ze(20),ax(20),zo(20),bx(20) common /block3/ h0,alfa,kswit common /block5/csph,snph,qcs,qsn,csps,snps common /block6/ ct,st c c change input to cylindrical coordinates rho = rr * st z = rr * ct c fl = 0.5 * (1.-alfa) c hhr = 0.0 hhz = 0.0 hhp = 0.0 thr = 0.0 thz = 0.0 d = 1.0d-02 c c do 100 i = 1,15 c arx = ze(i)/r * rho arz = ze(i)/r * z tarx = zo(i)/r * rho tarz = zo(i)/r * z c ex = dexp(+arz) ex = ex * fl * ax(i) tex = dexp(+tarz) tex = tex * fl * bx(i) c b0 = 1.0 b1 = 0.0 b2 = 0.0 tb0 = 1.0 tb1 = 0.0 if(abs(rho).ge.1.0d-06) then b0 = vbfj0(arx) b1 = vbfj1(arx) b2 = vbfj2(b1,b0,arx) tb0 = vbfj0(tarx) tb1 = vbfj1(tarx) end if hhr = hhr + (b0-b2)*ex hhp = hhp + (b0+b2)*ex hhz = hhz + 2.0*b1*ex c thr = thr + tb1*tex thz = thz + tb0*tex c c 100 continue c hhr = hhr*csps*csph hhp = hhp*csps* snph hhz = hhz*csps*csph thr = thr*snps thz = thz*snps c c hhr = - hhr + thr hhp = + hhp hhz = - hhz - thz c output in spherical coordinates hr = hhr*st + hhz*ct ht = hhr*ct - hhz*st hp = hhp c return end c --------------v btasph v---------------------------------------------- c --------------v vatan3 v---------------------------------------------- real*8 function vatan3(y,x) implicit double precision(a-h,o-z) c pi = 3.1415926535 vatan3 = 0.00 if((x.ne.0.0).or.(y.ne.0.0)) vatan3 = atan2(y,x) c return end c -------------v end of vatan3 v---------------------------------------- c -------------v vbfj0 v------------------------------------------------ real*8 function vbfj0(x) implicit double precision(a-h,o-z) c c purpose- calculate bessel function j0(x) c using 9.4.1 and 9.4.3 of ref2 to order 12 approximation c input- x, argument of the bessel function c output- vbfj0(x), =j0(x) c data a1/-2.2499997/,a2/1.2656208/,a3/-0.3163866/,a4/0.0444479/, 1 a5/-0.0039444/,a6/0.00021/,b0/0.79788456/,b1/-0.00000077/, 2 b2/-0.0055274/,b3/-0.00009512/,b4/0.00137237/,b5/-0.00072805/, 3 b6/0.00014476/,c0/-0.78539816/,c1/-0.04166397/,c2/-0.00003954/, 4 c3/0.00262573/,c4/-0.00054125/,c5/-0.00029333/,c6/0.00013558/ c x0=abs(x) t=x0/3.0 c if (t.lt.1.0) then c use eq. 9.4.1 t=t*t h1=1.0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6))))) vbfj0=h1 else c use eq. 9.4.3 t=1.0/t h2 =b0+t*(b1+t*(b2+t*(b3+t*(b4+t*(b5+t*b6))))) h3=x0+c0+t*(c1+t*(c2+t*(c3+t*(c4+t*(c5+t*c6))))) vbfj0=h2*dcos(h3)/dsqrt(x0) end if return end c -------------v end of vbfj0 v----------------------------------------- c -------------v vbfj1 v------------------------------------------------ real*8 function vbfj1(x) implicit double precision(a-h,o-z) c c purpose- calculate bessel function j1(x) c using 9.4.4 and 9.4.6 ref2 to order 12 approximation c input- x, argument of bessel function c output- bfj1, bessel function j1(x) c data a1/-.56249985/, 1 a2/.21093573/,a3/-.03954289/,a4/.00443319/,a5/-.00031761/, 2 a6/.00001109/,b0/.79788456/,b1/.00000156/,b2/.01659667/, 3 b3/.00017105/,b4/-.00249511/,b5/.00113653/,b6/-.00020033/, 4 c0/-2.35619449/,c1/.12499612/,c2/.0000565/,c3/-.00637879/, 5 c4/.00074348/,c5/.00079824/,c6/-.00029166/ c x0=dabs(x) x1=dsign(1.0d0,x) t=x0/3.0 c if (t.lt.1.0) then c use 9.4.4 t=t*t h1=0.5+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*a6))))) vbfj1=x*h1 else c use 9.4.6 t=1.0/t h2=b0+t*(b1+t*(b2+t*(b3+t*(b4+t*(b5+t*b6))))) h3=x0+c0+t*(c1+t*(c2+t*(c3+t*(c4+t*(c5+t*c6))))) vbfj1=x1*h2*dcos(h3)/dsqrt(x0) end if return end c -------------v end of vbfj1 v----------------------------------------- c -------------v vbfj2 v------------------------------------------------ real*8 function vbfj2(b1,b0,x) implicit double precision(a-h,o-z) c c purpose- calculate bessel function j2(x) using 9.1.27 ref2 c input- b1, j1(x) bessel from function bfj1 c b0, j0(x) bessel from function bfj0 c x, argument of bessel functions c output- bfj2, bessel j2(x) c vbfj2 = 0.0 if(x.gt.0.0) then vbfj2 = (2./x) * b1 - b0 end if c return end c -------------v end of vbfj2 v----------------------------------------- c -------------v cfita v------------------------------------------------ subroutine cfita(rr,hr,hz,hp) implicit double precision(a-h,o-z) c c purpose- calculate chapman-ferraro field of imf outside cylinder c input- rr rho in cylindrical coordinates c output- hr,hz,hp cylindrical field components c variables- himf, strength of imf field c xi, imf angle c r, radius of magnetopause sphere and tail c common /block1/ bex, r, q, pi, psi common /fimf/ himf, cdipol, cimf, xi c c get imf field components call bimf(himf,xi, d1,d2,d3,d4,d5,d6,d7, hr3,hz3,hp3) c fr = r/rr fr = fr*fr*fr c c components derived from 5.30 ref1 hr = - hr3 * fr hz = + 0.5 * hz3 * fr hp = + 0.5 * hp3 * fr c return end c -------------v end of cfita v------------------------------------------ c --------------v fun2 v------------------------------------------------- subroutine fun2(n,ndim) implicit double precision(a-h,o-z) c c purpose- calculate integrand for fourier-bessel integrals. one c for each of ndim grid points. done for tilt 0 and 90 deg c input- n, order of bessel function c ndim, number of grid points c output- y1,y0 integrands for tilt 0, 90 degrees, 5.24,5.25 ref1 c output to block 9 c variables- ze, zeros of derivatives of bessel j1 c zo, " " " " " j0 c x, radial values of grid points c r, radius of tail cylinder c common /block8/ w1(1001),w0(1001) common /block9/ x(1001),y1(1001),y0(1001),z1(1001),z0(1001), hs common /block1/ b ,r, q, pi, psi common /block2/ ze(20), ax(20), zo(20), bx(20) c c do 1 i = 1, ndim c arguments of bessels arx = ze(n) * x(i) / r tarx = zo(n) * x(i) / r c set bessel values for very small arguments b1 = 0.0 b0 = 1.0 if(abs(x(i)).ge.1.0d-06) then c get bessel function values j1 and j0 b1 = vbfj1(arx) b0 = vbfj0(tarx) end if c integrands of 5.24 and 5.25 ref1 y1(i) = x(i) * w1(i) * b1 y0(i) = x(i) * w0(i) * b0 1 continue c c return end c -------------v end of fun2 v------------------------------------- c -------------v vlegen v------------------------------------------ subroutine vlegen(t,ne) implicit double precision(a-h,o-z) c c purpose- calculate legendre polynomial terms c input- t, theta of spherical system c ne, order of legendre polynomial c output- pn, legendre polynomial terms m=0 c cpn, derivative of pn w.r.t. dcos(t) c pnm, legendre polynomial terms m=1 c common /vleg1/ pn(20), pnm(20), cpn(20) common /block6/ c1, s1 c c1 = dcos(t) s1 = dsin(t) c c2 = c1*c1 - s1*s1 s2 = 2.*s1*c1 c pn(1) = c1 pn(2) = 0.25*(3.*c2+1.) c cpn(1) = 1.0 cpn(2) = 3.*c1 c pnm(1) = + s1 * cpn(1) pnm(2) = + s1 * cpn(2) c c nele = ne if(ne.lt.3) return if(ne.eq.3) nele= 4 c if(ne.gt.3) nele= ne c do 1 n=3,nele fn = float(n) zfe = 2.*fn - 1.0 c recurrence relation from ref3 Jahnke-emde normalized pn(n) = (1./fn) * (zfe*c1*pn(n-1) - (fn-1.)*pn(n-2)) cpn(n) = zfe*pn(n-1) + cpn(n-2) pnm(n) = + s1*cpn(n) 1 continue c c return end c ------------v end of vlegen v---------------------------------------- c -------------v bata v------------------------------------------------ subroutine btata(z,rho,hhr,hhz,hhp) implicit double precision(a-h,o-z) c c purpose- calculate components of tailfield solution for zv>0 c iincluding stretch,summaries 6.3 and 6.4 ref1 c input- z,rho cylindrical coordinates c output- hhr,hhz,hhp cylindrical field components c variables- alfa, plasma content stretch parameter c snph,csph sin and cos of phi, a cylindrical coordinate c snps,csps sin and cos of psi, the dipole tilt angle c ax,bx coefficients for fourier-bessel series c r, radius of tail c ze,zo zeros of bessel function derivatives 5.25 ref1 c common /block1/ b,r,q,pi,psi common /block2/ ze(20),ax(20),zo(20),bx(20) common /block5/csph,snph,qcs,qsn,csps,snps common /block6/ ct,st common /block3/ h0,alfa,kswit c hhr = 0.0 hhz = 0.0 hhp = 0.0 thr = 0.0 thz = 0.0 d = 1.0d-02 c if(abs(alfa).lt.1.0d-06) then alfa = 1.0d-06 end if c fl = 0.5 * (1.-alfa) c c do 100 i = 1,15 c arguments of bessel functions arx = ze(i)/r*rho tarx = zo(i)/r*rho c arz = ze(i)/r*z tarz = zo(i)/r*z c flzx = alfa * arz flzy = alfa * tarz c ex = 0.0 tex = 0.0 if((arz.lt.170.).and.(tarz.lt.170.)) then ex = dexp(-arz) tex = dexp(-tarz) c ex = ex * fl * ax(i) tex = tex * fl * bx(i) end if c flex = 0.0 fley = 0.0 c filter changed from 170. to 160. if ((flzx.lt.160.).and.(flzy.lt.160.)) then flex = dexp(-flzx) fley = dexp(-flzy) end if c flex = flex * alfa * ax(i) fley = fley * alfa * bx(i) c b0 = 1.0 b1 = 0.0 b2 = 0.0 tb0 = 1.0 tb1 = 0.0 c if(abs(rho).ge.1.0d-06) then c get bessel function b0 = vbfj0(arx) b1 = vbfj1(arx) b2 = vbfj2(b1,b0,arx) tb0 = vbfj0(tarx) tb1 = vbfj1(tarx) end if c hhr = hhr + (b0-b2) * (flex + ex) hhp = hhp + (b0+b2) * (flex + ex) hhz = hhz + 2.*b1 * (flex/alfa + ex) c thr = thr + tb1 * (fley + tex) thz = thz + tb0 * (fley/alfa + tex) c 100 continue c superposition of 0 and 90 degree tilt cases hhr = hhr*csps*csph hhp = hhp*csps* snph hhz = hhz*csps*csph thr = thr*snps thz = thz*snps c c hhr = +hhr -thr hhp = -hhp hhz = -hhz -thz c return end c -----------v end of btata v------------------------------------------ c -----------v vacta v------------------------------------------------- subroutine vacta(z,rr,hhr,hhz,hhp) implicit double precision(a-h,o-z) c c purpose- calculate nightside (zv>0) vacuum field components of c tail solution from potential 5.26 ref1. solution c includes dipole plus chapman-ferraro. first called on c second iteration as first iter. has analytic coefficient c input- z,rr z and rho of cylindrical system c output- hhr,hhz,hhp field components in cylindrical system c variables- ze,zo zeros of derivatives of bessels j1 and j0 c ax,bx coefficients of fourier-bessel series from tailco c 5.24 and 5.25 ref1 c r, radius of tail c snps,csps sin and cos of psi, the dipole tilt angle c snph,snph sin and cos of phi, a cylindrical coordinate c c common /block1/ b,r,q,pi,psi common /block2/ ze(20),ax(20),zo(20),bx(20) common /block5/csph,snph,qcs,qsn,csps,snps common /block6/ ct,st c hhr = 0.0 hhz = 0.0 hhp = 0.0 thr = 0.0 thz = 0.0 d = 1.0d-02 c c do 100 i = 1,15 c arguments for bessel functions arx = ze(i)/r*rr arz = ze(i)/r*z tarx = zo(i)/r*rr tarz = zo(i)/r*z c ex = 0.0 tex = 0.0 c if((arz.lt.170.).and.(tarz.lt.170.)) then ex = dexp(-arz) ex = ex*ax(i) tex = dexp(-tarz) tex = tex*bx(i) end if c b0 = 1.0 b1 = 0.0 b2 = 0.0 tb0 = 1.0 tb1 = 0.0 c if(abs(rr).ge.1.0d-06) then b0 = vbfj0(arx) b1 = vbfj1(arx) b2 = vbfj2(b1,b0,arx) tb0 = vbfj0(tarx) tb1 = vbfj1(tarx) end if c hhr = hhr + (b0-b2)*ex hhp = hhp + (b0+b2)*ex hhz = hhz + 2.0*b1*ex c thr = thr + tb1*tex thz = thz + tb0*tex c 100 continue c superimpose 0 and 90 degree tilt cases hhr = hhr*csps*csph hhp = hhp*csps* snph hhz = hhz*csps*csph thr = thr*snps thz = thz*snps c c hhr = +hhr -thr hhp = -hhp hhz = -hhz -thz c return end c ------------v end of vacta v------------------------------------------ c ------------v tapot v------------------------------------------------- subroutine tapot(ndim) implicit double precision(a-h,o-z) c c purpose- get sum of chapman-ferraro and dipole potentials at ndim c locations on the cylinder cap zv=0 c input- ndim, number of grid points for integration c output- w1,wo sum of potentials for tilts 0 and 90 degree c output is potential*(1/dipole moment0 c variables- hs, grid size c r, radius of tail c x, radial positions of grid points c common /block8/ w1(1001),w0(1001) common /block9/ x(1001),y1(1001),y0(1001),z1(1001),z0(1001), hs common /block1/ b, r, q, pi, psi c c generate grid points on cylinder cap x(1) = 0.0 hs = r/float(ndim-1) do 1 i = 2, ndim x(i) = float(i-1) * hs 1 continue c c get potentials of chapman-ferraro tilt 0 and 90 degrees c eq. 5.18-5.20 ref1 do 2 i = 1, ndim call cfpot(x(i),pot1,pot0) w1(i) = pot1 w0(i) = pot0 2 continue c c get potentials of dipole,tilt 0 and 90 degrees 5.17 ref1 do 3 i = 1, ndim call dipot(x(i),pot1,pot0) w1(i) = w1(i) + pot1 w0(i) = w0(i) + pot0 3 continue c c return end c -------------v end of tapot v---------------------------------------- c -------------v cfpot v----------------------------------------------- subroutine cfpot(rr,pot1,pot0) implicit double precision(a-h,o-z) c c purpose- calculate chapman-ferraro potential for a point at the c cylinder lid, zv=0. uses 5.11 with theta=90 degrees, c equivalent to 5.18-5.20 ref1. both tilt cases 0 and 90 c degrees done. c input- rr, rho in cylindrical, distance from tail axis c output- pot1, potential for 0 degree tilt case c pot0, potential for 90 degree tilt case c output is actually potential*(-1/dipole moment) c variables- qr, dipole moment c p1,po legendre polynomial terms for m=1 and 0, respect. c g1,go coefficients of legendre terms, m=1 and 0 respect. c r, radius of tail c common /block1/ b, r, q, pi, psi common /chablo/ g1(20), g0(20) c dimension p0(21), p1(21) c c legendre terms with theta=90 degrees p0(2) = - 0.5 p1(1) = + 1.0 p1(3) = 3.*p0(2) c c do 1 ny = 2, 10 c every other term is zero, recurrence relation used c subrout. legen not called zfny = 2.*float(ny) fakny = (zfny-1.)/zfny p0(2*ny) = - p0(2*ny-2) * fakny p1(2*ny+1) = + (zfny+1.) * p0(2*ny) 1 continue c c tilt = 0 degrees case rrp = rr rp = (rr*rr)/(r*r) c pot1 = g1(1) * p1(1) * rr c do 2 ny = 1, 9 zfe = float(2*ny+1) rrp = rrp * rp c safety switch on 7/1/93 ............................. if(rrp.le.1.0d-30) goto 4 if(dlog(rrp).lt.-25.0*dlog(10.0d0)) goto 2 pot1 = pot1 + g1(2*ny+1)*p1(2*ny+1)*rrp/zfe 2 continue 4 continue c c do 2 ny = 1, 9 c zfe = float(2*ny+1) c rrp = rrp * rp c pot1 = pot1 + g1(2*ny+1) * p1(2*ny+1) * rrp / zfe c 2 continue c qr = q * (r**3) pot1 = -pot1/qr c c tilt = 90 degree case rrp = r rp = (rr*rr)/(r*r) c pot0 = 0.0 c do 3 ny = 1, 10 zfe = float(2*ny) rrp = rrp * rp c safety switch added on 7/1/93 ....................... if(rrp.le.1.0d-30) goto 5 if(dlog(rrp).lt.-30.0*dlog(10.0d0))go to 3 pot0 = pot0 + g0(2*ny) * p0(2*ny) * rrp / zfe 3 continue 5 continue c do 3 ny = 1, 10 c zfe = float(2*ny) c rrp = rrp * rp c pot0 = pot0 + g0(2*ny) * p0(2*ny) * rrp / zfe c 3 continue c qr = q * (r**3) pot0 = -pot0/qr c return end c -------------v end of cfpot v----------------------------------------- c -------------v dipot v------------------------------------------------ subroutine dipot(rr,pot1,pot2) implicit double precision(a-h,o-z) c c purpose- calculate dipole potential at lid of cylinder using 5.17 c ref1 for tilt 0 and 90 degree cases c input- rr, rho in cylindrical, distance from tail axis c b, = radius minus stand-off distance c output- pot1,pot2 dipole potential for tilt 0 and 90 degrees c output is potential*(-1/dipole moment) c common /block1/ b, r, q, pi, psi c rrb = dsqrt(rr*rr + b*b) rrb = rrb**3 c pot1 = + rr / rrb pot2 = + b / rrb c return end c --------------v end of dipot v------------------------------------- c --------------v trans1 v------------------------------------------- subroutine trans1(grr,gthet,gphi,xd,yd,zd) implicit double precision(a-h,o-z) rr = grr fac = 3.141592654/180.0 thet = gthet*fac pphi = gphi*fac sst = dsin(thet) cct = dcos(thet) ssp = dsin(pphi) ccp = dcos(pphi) xd = rr*sst*ccp yd = rr*sst*ssp zd = rr*cct return end c -------------v end of trans1 v----------------------------------------- c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c !!!!! The End !!!!! c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!