| 1 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2 | subroutine sol2ls(sol,ls) |
|---|
| 3 | |
|---|
| 4 | c convert a given martian day number (sol) |
|---|
| 5 | c into corresponding solar longitude, Ls (in degr.), |
|---|
| 6 | c where sol=0=Ls=0 is the |
|---|
| 7 | c northern hemisphere spring equinox. |
|---|
| 8 | |
|---|
| 9 | implicit none |
|---|
| 10 | |
|---|
| 11 | c input |
|---|
| 12 | real sol |
|---|
| 13 | c output |
|---|
| 14 | real ls |
|---|
| 15 | |
|---|
| 16 | c local variables |
|---|
| 17 | double precision year_day,peri_day,timeperi,e_elips |
|---|
| 18 | double precision pi,radtodeg |
|---|
| 19 | c number of martian days (sols) in a martian year |
|---|
| 20 | parameter (year_day=668.6d0) |
|---|
| 21 | c perihelion date (in sols) |
|---|
| 22 | parameter (peri_day=485.35d0) |
|---|
| 23 | c orbital eccentricity |
|---|
| 24 | parameter (e_elips=0.09340d0) |
|---|
| 25 | parameter (pi=3.14159265358979d0) |
|---|
| 26 | c radtodeg: 180/pi |
|---|
| 27 | parameter (radtodeg=57.2957795130823d0) |
|---|
| 28 | c timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99 |
|---|
| 29 | parameter (timeperi=1.90258341759902d0) |
|---|
| 30 | |
|---|
| 31 | double precision zanom,xref,zx0,zdx,zteta,zz |
|---|
| 32 | c xref: mean anomaly, zx0: eccentric anomaly, zteta: true anomaly |
|---|
| 33 | integer iter |
|---|
| 34 | |
|---|
| 35 | zz=(sol-peri_day)/year_day |
|---|
| 36 | zanom=2.*pi*(zz-nint(zz)) |
|---|
| 37 | xref=dabs(zanom) |
|---|
| 38 | |
|---|
| 39 | c The equation zx0 - e * sin (zx0) = xref, solved by Newton |
|---|
| 40 | zx0=xref+e_elips*dsin(xref) |
|---|
| 41 | do 110 iter=1,10 |
|---|
| 42 | zdx=-(zx0-e_elips*dsin(zx0)-xref)/(1.-e_elips*dcos(zx0)) |
|---|
| 43 | if(dabs(zdx).le.(1.d-7)) then ! typically, 2 or 3 iterations are enough |
|---|
| 44 | goto 120 |
|---|
| 45 | endif |
|---|
| 46 | zx0=zx0+zdx |
|---|
| 47 | 110 continue |
|---|
| 48 | 120 continue |
|---|
| 49 | zx0=zx0+zdx |
|---|
| 50 | if(zanom.lt.0.) zx0=-zx0 |
|---|
| 51 | |
|---|
| 52 | c compute true anomaly zteta, now that eccentric anomaly zx0 is known |
|---|
| 53 | zteta=2.*datan(dsqrt((1.+e_elips)/(1.-e_elips))*dtan(zx0/2.)) |
|---|
| 54 | |
|---|
| 55 | c compute Ls |
|---|
| 56 | ls=real(zteta-timeperi) |
|---|
| 57 | if(ls.lt.0.) ls=ls+2.*real(pi) |
|---|
| 58 | if(ls.gt.2.*pi) ls=ls-2.*real(pi) |
|---|
| 59 | c convert Ls in deg. |
|---|
| 60 | ls=real(radtodeg)*ls |
|---|
| 61 | |
|---|
| 62 | return |
|---|
| 63 | end |
|---|
| 64 | |
|---|
| 65 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 66 | real function ls2sol(ls) |
|---|
| 67 | |
|---|
| 68 | c Returns solar longitude, Ls (in deg.), from day number (in sol), |
|---|
| 69 | c where sol=0=Ls=0 at the northern hemisphere spring equinox |
|---|
| 70 | |
|---|
| 71 | implicit none |
|---|
| 72 | |
|---|
| 73 | c Arguments: |
|---|
| 74 | real ls |
|---|
| 75 | |
|---|
| 76 | c Local: |
|---|
| 77 | double precision xref,zx0,zteta,zz |
|---|
| 78 | c xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly |
|---|
| 79 | double precision year_day |
|---|
| 80 | double precision peri_day,timeperi,e_elips |
|---|
| 81 | double precision pi,degrad |
|---|
| 82 | parameter (year_day=668.6d0) ! number of sols in a amartian year |
|---|
| 83 | c data peri_day /485.0/ |
|---|
| 84 | parameter (peri_day=485.35d0) ! date (in sols) of perihelion |
|---|
| 85 | c timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99 |
|---|
| 86 | parameter (timeperi=1.90258341759902d0) |
|---|
| 87 | parameter (e_elips=0.0934d0) ! eccentricity of orbit |
|---|
| 88 | parameter (pi=3.14159265358979d0) |
|---|
| 89 | parameter (degrad=57.2957795130823d0) |
|---|
| 90 | |
|---|
| 91 | if (abs(ls).lt.1.0e-5) then |
|---|
| 92 | if (ls.ge.0.0) then |
|---|
| 93 | ls2sol = 0.0 |
|---|
| 94 | else |
|---|
| 95 | ls2sol = real(year_day) |
|---|
| 96 | end if |
|---|
| 97 | return |
|---|
| 98 | end if |
|---|
| 99 | |
|---|
| 100 | zteta = ls/degrad + timeperi |
|---|
| 101 | zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips))) |
|---|
| 102 | xref = zx0-e_elips*dsin(zx0) |
|---|
| 103 | zz = xref/(2.*pi) |
|---|
| 104 | ls2sol = real(zz*year_day + peri_day) |
|---|
| 105 | if (ls2sol.lt.0.0) ls2sol = ls2sol + real(year_day) |
|---|
| 106 | if (ls2sol.ge.year_day) ls2sol = ls2sol - real(year_day) |
|---|
| 107 | |
|---|
| 108 | return |
|---|
| 109 | end |
|---|
