C****************************************************************************** SUBROUTINE csheet_struc(ZNS3,R,theta,phi,XJSO,YJSO,ctime) C !Version 2004 08 15 C INPUT: R,theta,phi ->in System III C XJSO,YJSO ->in JSO C ctime ->UNIX cline time C OUTPUT: ZNS3 ->height of the current sheet in System III C IMPLICIT REAL*8(A-H,O-Z) REAL*8, INTENT(IN) :: R,theta,phi,XJSO,YJSO,ctime REAL*8, INTENT(OUT) :: ZNS3 DATA XH /47/ DATA dLamVIP4 /5.92068/ DATA dThetaVIP4 /0.16615534478986/ RHO=R*dsin(theta) CALL JSun(ctime,stheta,sphi,phase) RLT=datan2(YJSO,XJSO) !zero at noon dlamprime=waveDelay(RLT,RHO)+bendBack(RLT,RHO)+dLamVIP4 rho1=dsqrt((XH*dtanh(dabs(XJSO/XH)))**2+YJSO**2) term2=rho1*dtan(dThetaVIP4) ZNS3=term2*dcos(phi-dlamprime)+XJSO*(1.0-dtanh(dabs(XH/XJSO))) ~ *dtan(stheta) RETURN END C****************************************************************************** C****************************************************************************** REAL*8 FUNCTION waveDelay(RLT,RHO) IMPLICIT REAL*8(A-H,O-Z) REAL*8, INTENT(IN) :: RLT, RHO Data B0 / 3.79406064748764D-002/ Data C1 / 1.34045565128326 / !C1=Domegaj*rho2/V0 !with Domegaj=36.26125*pi/180 Data rho2/ 83.3625411987305/ Data B1/ 0.748919927647621/, Phi1/ -1.04848085782547/ Data B2/ 0.614373903493367/, Phi2/ -1.95757783830202/ Data B3/ 0.541435668088830/, Phi3/ -1.25514203423947/ Data B4/ 7.02600475148873D-002/, Phi4/ -1.25542899066018/ Data B5/ 26.1889910697937/ rhoor2=RHO/rho2 rhooB5=RHO/B5 waveDelay=B0-C1*dlog(dcosh(rhoor2))* ~(1+dtanh(rhooB5)*(B1*dcos(RLT-Phi1)+B2*dcos(2*RLT-Phi2)+ ~B3*dcos(3*RLT-Phi3)+B4*dcos(4*RLT-Phi4))) RETURN END C****************************************************************************** C****************************************************************************** REAL*8 FUNCTION bendBack(RLT,RHO) IMPLICIT REAL*8(A-H,O-Z) REAL*8 RLT, RHO Data rho0 /33.0/ Data A1 / 1.062528141129624e-004 /!divided by rho1 (1.6*rho0=52.8) Data A2 / 5.950493765539412e-005 /!divided by rho1 (1.6*rho0=52.8) Data A3 / 6.788463593277820e-005 /!divided by rho1 (1.6*rho0=52.8) Data A4 / 3.009364679056859e-005 /!divided by rho1 (1.6*rho0=52.8) Data Phi1 / 1.79696364159234/ Data Phi2 / 0.41408648226213/ Data Phi3 /-0.22365844912454/ Data Phi4 /-0.42846193975833/ Data Cb1n /-0.18591516533271 /!Cb1=Cb1*rho0 Data Cb1logc/-8.534743746501581e-005/!Cb1n*dlog(dcosh(dabs(1/rho0))) Data ddK2 /-0.01893939393939 /!dK2=-1/(1.6*rho0) ddK1=A1*dcos( RLT-Phi1) ddK3=A2*dcos(2*RLT-Phi2) ddK5=A3*dcos(3*RLT-Phi3) ddK7=A4*dcos(4*RLT-Phi4) T3=dexp(ddK2*rho)/(ddK2**2)*(ddK2*rho-1) - ~ dexp(ddK2 )/(ddK2**2)*(ddK2 -1) T4=ddK1/2+ddK3/2+ddK5/2+ddK7/2 rhoSq=rho**2 bendBack=Cb1n*dlog(dcosh(dabs(rho/rho0))) ~ -Cb1logc+T4*rhoSq+T3*(-ddK1-ddK3-ddK5-ddK7)-T4 RETURN END C****************************************************************************** c****************************************************************************** Subroutine JSun(ctime,stheta,sphi,phase) c INPUT: ctime of the data point C OUTPUTS: stheta, sphi, latitude and longitude (in radians) of the Sun in system III (RH). C OUTPUTS: phase, Orbital phase angle of Jupiter (in radians) c The equations are written in etime the J200 time convention followed by PDS. c We first convert ctime to etime c Last updated August 12, 2002. c C The program first calculates the direction of the sun in non-rotating C coordinates.from equations of the type: c theta=a1*cos(omegay*t)+a2*sin(omegay*t)+a3*cos(2.*omegay*t)+ c a4*sin(2.*omegay*t)+a5 c fphi=b1*cos(omegay*t)+b2*sin(omegay*t)+b3*cos(2.*omegay*t)+ c b4*sin(2.*omegay*t)+b5 c Then we rotate into the System III coordinates c ctime is in double precision. theta and phi are single precision c +variables. c c c omega is jupiter's rotation rate c omegay is jupiter's yearly orbital rate c Initialize variables IMPLICIT REAL*8(A-H,O-Z) Dimension aa(7), bb(7), x(7) REAL*8, INTENT(IN) :: ctime REAL*8 :: etimer,e_time1,fphi,yrjup,omega,D360,omegay,t,year REAL*8, INTENT(OUT) :: stheta,sphi,phase Parameter (PI=3.1415927,twopi=2.*PI,radian=PI/180.,degree=180./PI) Parameter (yrjup=.1185652502D2*86400d0*.36525d3) Parameter (omega=870.536/86400.0,omegay=2.*PI/yrjup) Parameter (year=86400D0*.36525D3,D360=360.0) Parameter (e_time1=-8.25767955817479d8) Parameter (three=3.123*radian) Parameter (tan3=0.054560676,sin3=0.054479647,cos3=0.99851488) Data aa /0.14347029,3.1145815,-0.12025561,0.093909436, ~-0.39321884e-5,0.10194945e-3,-0.12799464/ Data bb /-4.5467523,3.1848875,-0.16329986,-0.09776818, ~0.17556527e-3,-0.01978317,44.55915/ c First calculate the latitude and longitude in non-rotating Jupiter coordinates. c Calculate the best fit theta and fphi t = etimer(ctime) - e_time1 x(1)=dcos(omegay*t) x(2)=dsin(omegay*t) x(3)=dcos(2.*omegay*t) x(4)=dsin(2.*omegay*t) x(5)=(t/year)**2 x(6)=t/year x(7)=1.0 stheta=0. c fphi is phi in Jupiter fixed (non-rotating) coordinate fphi=0.0 Do 1 j=1,7 fphi=fphi+bb(j)*x(j) stheta=stheta+aa(j)*x(j) 1 continue C Now rotate the longitude to Jupiter System III c First Add the rotation of Jupiter around the Sun. c fphi is the phi of the Sun as unspinning Jupiter goes around the sun fphi=DMod(fphi+t/yrjup*360.d0, D360) c Next add the rotation of Jupiter around its axis. sphi=DMod(fphi-t*omega, D360) if (sphi .lt. 0.0) sphi = sphi+360.0 sphi=sphi*radian stheta=stheta*radian c Now compute the orbital phase (called phi2 or phase here) c There are two solutions to the problem. Only one that is close to phi2b is correct. if(stheta .ge. three)stheta=(three) if(-stheta .ge. three)stheta=-(three) phi21=(sphi+PI)+dacos(dtan(stheta)/tan3) phi22=(sphi+PI)-dacos(dtan(stheta)/tan3) phi21=dMOD(phi21,twopi) phi22=dMOD(phi22,twopi) if(phi22 .lt. 0) phi22=phi22+twopi dphi=fphi+48.23012 phi2b=sphi-dphi*radian phi2b=dMOD(phi2b,twopi) if(phi2b .lt. 0) phi2b=phi2b+twopi phase=phi21 dif2=dabs(phi22-phi2b) if(dif2 .gt. 350*radian) dif2=twopi-dif2 dif1=dabs(phi21-phi2b) if(dif1 .gt. 350*radian) dif1=twopi-dif1 if(dif2 .lt. dif1) phase=phi22 c phase=phi2b phase=dMOD(phase,twopi) c Write(*,*)stheta*degree,phase*degree,phi2b*degree Return end c******************************************************************* c******************************************************************* REAL*8 Function etimer(ctime) Real*8, INTENT(IN) :: ctime REAL*8 Tcor C Integer ctime, etime,Tcor Tcor=0 if(ctime .ge. 189302400) Tcor = Tcor+10 if(ctime .ge. 205027200) Tcor = Tcor+1 if(ctime .ge. 220924800) Tcor = Tcor+1 if(ctime .ge. 252460800) Tcor = Tcor+1 if(ctime .ge. 283996800) Tcor = Tcor+1 if(ctime .ge. 315532800) Tcor = Tcor+1 if(ctime .ge. 347155200) Tcor = Tcor+1 if(ctime .ge. 378691200) Tcor = Tcor+1 if(ctime .ge. 410227200) Tcor = Tcor+1 if(ctime .ge. 441763200) Tcor = Tcor+1 if(ctime .ge. 489024000) Tcor = Tcor+1 if(ctime .ge. 520560000) Tcor = Tcor+1 if(ctime .ge. 552096000) Tcor = Tcor+1 if(ctime .ge. 615254400) Tcor = Tcor+1 if(ctime .ge. 694224000) Tcor = Tcor+1 if(ctime .ge. 757382400) Tcor = Tcor+1 if(ctime .ge. 788918400) Tcor = Tcor+1 if(ctime .ge. 836179200) Tcor = Tcor+1 if(ctime .ge. 867715200) Tcor = Tcor+1 if(ctime .ge. 899251200) Tcor = Tcor+1 if(ctime .ge. 946684800) Tcor = Tcor+1 if(ctime .ge. 993945600) Tcor = Tcor+1 if(ctime .ge. 1041379200) Tcor = Tcor+1 etimer=ctime + dble(Tcor) - .1072958367816D10 C etime = ctime + Tcor - 1072958367 Return END C******************************************************************** Real*8 Function ctimer(iyr,imon,iday,ihr,imin,sec) integer, INTENT(IN) :: iyr,imon,iday,ihr,imin integer*4 ndays,doy Real*8, INTENT(IN) :: sec c First calculate the number of days from Jan 1, 1966 ndays=0 If (iyr .ge. 1966) then Do 1 i=1966,iyr-1 ndays=ndays+365 if(((mod(i,4) .eq. 0) .and. (mod(i,100) .ne. 0)) .or. + (mod(i,400)) .eq. 0) ndays=ndays+1 1 continue c Now add the number of days of the current year ndays=ndays+doy(iyr,imon,iday)-1 ctimer = dble(ndays)*86400D0+dble(ihr)*3600D0+dble(imin)*60D0+ + dble(sec) go to 4 end if c Calculate the seconds for the negative years If (iyr .lt. 1966) then Do 2 i=iyr, 1965 ndays=ndays-365 if(((mod(i,4) .eq. 0) .and. (mod(i,100) .ne. 0)) .or. + (mod(i,400)) .eq. 0) ndays=ndays-1 2 continue c Now subtract the number of days of the current year ndays=ndays+doy(iyr,imon,iday) end if ctimer = dble(ndays-1)*86400D0+dble(ihr)*3600D0+dble(imin)*60D0+ + dble(sec) 4 continue idyr=doy(iyr,imon,iday) Write(*,*) " doy is ",idyr Write(*,*)" ndays is ",ndays Return End c******************************************************************* Integer Function doy(iyr,imon,iday) Integer, INTENT(IN) :: iyr,imon,iday Integer mon(12) data mon /31,28,31,30,31,30,31,31,30,31,30,31/ doy=0 Do 1 i = 2,imon doy=doy+mon(i-1) c Add an extra day for February if(i .eq. 3) then if(((mod(iyr,4) .eq. 0) .and. (mod(iyr,100) .ne. 0)) .or. + (mod(iyr,400)) .eq. 0) doy=doy+1 end if 1 continue doy=doy+iday Return End c*******************************************************************