| 1 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 2 | real function ls2sol(ls) |
|---|
| 3 | |
|---|
| 4 | c Returns solar longitude, Ls (in deg.), from day number (in sol), |
|---|
| 5 | c where sol=0=Ls=0 at the northern hemisphere spring equinox |
|---|
| 6 | |
|---|
| 7 | implicit none |
|---|
| 8 | |
|---|
| 9 | c Arguments: |
|---|
| 10 | real ls |
|---|
| 11 | |
|---|
| 12 | c Local: |
|---|
| 13 | double precision xref,zx0,zteta,zz |
|---|
| 14 | c xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly |
|---|
| 15 | double precision year_day |
|---|
| 16 | double precision peri_day,timeperi,e_elips |
|---|
| 17 | double precision pi,degrad |
|---|
| 18 | parameter (year_day=668.6d0) ! number of sols in a amartian year |
|---|
| 19 | c data peri_day /485.0/ |
|---|
| 20 | parameter (peri_day=485.35d0) ! date (in sols) of perihelion |
|---|
| 21 | c timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99 |
|---|
| 22 | parameter (timeperi=1.90258341759902d0) |
|---|
| 23 | parameter (e_elips=0.0934d0) ! eccentricity of orbit |
|---|
| 24 | parameter (pi=3.14159265358979d0) |
|---|
| 25 | parameter (degrad=57.2957795130823d0) |
|---|
| 26 | |
|---|
| 27 | if (abs(ls).lt.1.0e-5) then |
|---|
| 28 | if (ls.ge.0.0) then |
|---|
| 29 | ls2sol = 0.0 |
|---|
| 30 | else |
|---|
| 31 | ls2sol = year_day |
|---|
| 32 | end if |
|---|
| 33 | return |
|---|
| 34 | end if |
|---|
| 35 | |
|---|
| 36 | zteta = ls/degrad + timeperi |
|---|
| 37 | zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips))) |
|---|
| 38 | xref = zx0-e_elips*dsin(zx0) |
|---|
| 39 | zz = xref/(2.*pi) |
|---|
| 40 | ls2sol = zz*year_day + peri_day |
|---|
| 41 | if (ls2sol.lt.0.0) ls2sol = ls2sol + year_day |
|---|
| 42 | if (ls2sol.ge.year_day) ls2sol = ls2sol - year_day |
|---|
| 43 | |
|---|
| 44 | return |
|---|
| 45 | end |
|---|
| 46 | |
|---|
| 47 | cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
