| 1 | MODULE solarlong_mod |
|---|
| 2 | |
|---|
| 3 | IMPLICIT NONE |
|---|
| 4 | |
|---|
| 5 | CONTAINS |
|---|
| 6 | |
|---|
| 7 | SUBROUTINE solarlong(pday,psollong) |
|---|
| 8 | use planete_h, only: lsperi, peri_day, year_day, e_elips |
|---|
| 9 | use comcstfi_h, only: pi |
|---|
| 10 | IMPLICIT NONE |
|---|
| 11 | |
|---|
| 12 | c======================================================================= |
|---|
| 13 | c Compute solar longitude psollong (in radians) for a given Mars date |
|---|
| 14 | c pday (where pday==0 at northern Spring Equinox and pday in sols |
|---|
| 15 | c and fractions thereof) |
|---|
| 16 | c======================================================================= |
|---|
| 17 | |
|---|
| 18 | c arguments: |
|---|
| 19 | c ---------- |
|---|
| 20 | |
|---|
| 21 | REAL,INTENT(IN) :: pday |
|---|
| 22 | REAL,INTENT(OUT) :: psollong |
|---|
| 23 | |
|---|
| 24 | c Local variables: |
|---|
| 25 | c ---------------- |
|---|
| 26 | |
|---|
| 27 | REAL :: xref ! mean anomaly |
|---|
| 28 | REAL :: zx0 ! eccentric anomaly |
|---|
| 29 | REAL :: zteta ! true anomaly |
|---|
| 30 | REAL :: zanom,zdx,zz |
|---|
| 31 | INTEGER :: iter |
|---|
| 32 | |
|---|
| 33 | |
|---|
| 34 | c compute mean anomaly |
|---|
| 35 | |
|---|
| 36 | zz=(pday-peri_day)/year_day |
|---|
| 37 | zanom=2.*pi*(zz-nint(zz)) |
|---|
| 38 | xref=abs(zanom) |
|---|
| 39 | |
|---|
| 40 | c solve equation zx0 - e * sin (zx0) = xref |
|---|
| 41 | c using Newton method |
|---|
| 42 | |
|---|
| 43 | zx0=xref+e_elips*sin(xref) |
|---|
| 44 | DO iter=1,10 ! 10 is overkill, typically converges in 1-3 iterations |
|---|
| 45 | zdx=-(zx0-e_elips*sin(zx0)-xref)/(1.-e_elips*cos(zx0)) |
|---|
| 46 | if(abs(zdx).le.(1.e-7)) exit |
|---|
| 47 | zx0=zx0+zdx |
|---|
| 48 | ENDDO |
|---|
| 49 | |
|---|
| 50 | zx0=zx0+zdx |
|---|
| 51 | if(zanom.lt.0.) zx0=-zx0 |
|---|
| 52 | |
|---|
| 53 | c compute true anomaly zteta, now that eccentric anomaly zx0 is known |
|---|
| 54 | |
|---|
| 55 | zteta=2.*atan(sqrt((1.+e_elips)/(1.-e_elips))*tan(zx0/2.)) |
|---|
| 56 | |
|---|
| 57 | c compute solar longitude from zteta |
|---|
| 58 | psollong=zteta+lsperi |
|---|
| 59 | |
|---|
| 60 | c handle limit cases where psollong lands outside [0:2*pi] |
|---|
| 61 | IF(psollong.LT.0.) psollong=psollong+2.*pi |
|---|
| 62 | IF(psollong.GT.2.*pi) psollong=psollong-2.*pi |
|---|
| 63 | |
|---|
| 64 | END SUBROUTINE solarlong |
|---|
| 65 | |
|---|
| 66 | END MODULE solarlong_mod |
|---|
