| 1 | def ls2sol(lstabin): |
|---|
| 2 | |
|---|
| 3 | # Returns solar longitude, Ls (in deg.), from day number (in sol), |
|---|
| 4 | # where sol=0=Ls=0 at the northern hemisphere spring equinox |
|---|
| 5 | import numpy as np |
|---|
| 6 | import math as m |
|---|
| 7 | |
|---|
| 8 | year_day=668.6 |
|---|
| 9 | peri_day=485.35 |
|---|
| 10 | timeperi=1.90258341759902 |
|---|
| 11 | e_elips=0.0934 |
|---|
| 12 | pi=3.14159265358979 |
|---|
| 13 | degrad=57.2957795130823 |
|---|
| 14 | |
|---|
| 15 | if type(lstabin).__name__ in ['int','float']: |
|---|
| 16 | lstab=[lstabin] |
|---|
| 17 | lsout=np.zeros([1]) |
|---|
| 18 | else: |
|---|
| 19 | lstab=lstabin |
|---|
| 20 | lsout=np.zeros([len(lstab)]) |
|---|
| 21 | |
|---|
| 22 | i=0 |
|---|
| 23 | for ls in lstab: |
|---|
| 24 | if (np.abs(ls) < 1.0e-5): |
|---|
| 25 | if (ls >= 0.0): |
|---|
| 26 | ls2sol = 0.0 |
|---|
| 27 | else: |
|---|
| 28 | ls2sol = year_day |
|---|
| 29 | else: |
|---|
| 30 | |
|---|
| 31 | zteta = ls/degrad + timeperi |
|---|
| 32 | zx0 = 2.0*m.atan(m.tan(0.5*zteta)/m.sqrt((1.+e_elips)/(1.-e_elips))) |
|---|
| 33 | xref = zx0-e_elips*m.sin(zx0) |
|---|
| 34 | zz = xref/(2.*pi) |
|---|
| 35 | ls2sol = zz*year_day + peri_day |
|---|
| 36 | if (ls2sol < 0.0): ls2sol = ls2sol + year_day |
|---|
| 37 | if (ls2sol >= year_day): ls2sol = ls2sol - year_day |
|---|
| 38 | |
|---|
| 39 | lsout[i]=ls2sol |
|---|
| 40 | i=i+1 |
|---|
| 41 | |
|---|
| 42 | return lsout |
|---|
| 43 | |
|---|
| 44 | def sol2ls(soltabin): |
|---|
| 45 | |
|---|
| 46 | # convert a given martian day number (sol) |
|---|
| 47 | # into corresponding solar longitude, Ls (in degr.), |
|---|
| 48 | # where sol=0=Ls=0 is the |
|---|
| 49 | # northern hemisphere spring equinox. |
|---|
| 50 | import numpy as np |
|---|
| 51 | import math as m |
|---|
| 52 | |
|---|
| 53 | year_day=668.6 |
|---|
| 54 | peri_day=485.35 |
|---|
| 55 | e_elips=0.09340 |
|---|
| 56 | radtodeg=57.2957795130823 |
|---|
| 57 | timeperi=1.90258341759902 |
|---|
| 58 | |
|---|
| 59 | if type(soltabin).__name__ in ['int','float']: |
|---|
| 60 | soltab=[soltabin] |
|---|
| 61 | solout=np.zeros([1]) |
|---|
| 62 | else: |
|---|
| 63 | soltab=soltabin |
|---|
| 64 | solout=np.zeros([len(soltab)]) |
|---|
| 65 | |
|---|
| 66 | i=0 |
|---|
| 67 | for sol in soltab: |
|---|
| 68 | zz=(sol-peri_day)/year_day |
|---|
| 69 | zanom=2.*np.pi*(zz-np.floor(zz)) |
|---|
| 70 | xref=np.abs(zanom) |
|---|
| 71 | |
|---|
| 72 | # The equation zx0 - e * sin (zx0) = xref, solved by Newton |
|---|
| 73 | zx0=xref+e_elips*m.sin(xref) |
|---|
| 74 | iter=0 |
|---|
| 75 | while iter <= 10: |
|---|
| 76 | iter=iter+1 |
|---|
| 77 | zdx=-(zx0-e_elips*m.sin(zx0)-xref)/(1.-e_elips*m.cos(zx0)) |
|---|
| 78 | if(np.abs(zdx) <= (1.e-7)): |
|---|
| 79 | continue |
|---|
| 80 | zx0=zx0+zdx |
|---|
| 81 | zx0=zx0+zdx |
|---|
| 82 | |
|---|
| 83 | if(zanom < 0.): zx0=-zx0 |
|---|
| 84 | # compute true anomaly zteta, now that eccentric anomaly zx0 is known |
|---|
| 85 | zteta=2.*m.atan(m.sqrt((1.+e_elips)/(1.-e_elips))*m.tan(zx0/2.)) |
|---|
| 86 | |
|---|
| 87 | # compute Ls |
|---|
| 88 | ls=zteta-timeperi |
|---|
| 89 | if(ls < 0.): ls=ls+2.*np.pi |
|---|
| 90 | if(ls > 2.*np.pi): ls=ls-2.*np.pi |
|---|
| 91 | # convert Ls in deg. |
|---|
| 92 | ls=radtodeg*ls |
|---|
| 93 | solout[i]=ls |
|---|
| 94 | i=i+1 |
|---|
| 95 | |
|---|
| 96 | return solout |
|---|
