! SPHER_HARM: Gravitational potential routines for satellites MODULE spher_harm USE fund_const USE output_control USE planet_masses IMPLICIT NONE PRIVATE ! ROUTINES ! geopot : Input of geopotential coefficients ! namcoe : Computes the l, m and i of an harmonic with a given ndex ! parm10 : Computes accelerations due to spherical harmonics ! (Not normalized) ! legit : Recursive computation of Legendre polynomials ! and Legendre associated functions ! ! FUNCTIONS ! indcoe : Address of harmonics (l,m,i) ! numcoe : Number of coefficients up to degree ites ! fnspar : Computes the factor for the normalization ! of spherical harmonic coefficients ! ! narmx=max.no.harmonic coefficients, itesx=max.degree and order INTEGER, PARAMETER, PUBLIC :: itesx=50 INTEGER, PARAMETER, PUBLIC :: narmx=(itesx+1)**2 INTEGER, PUBLIC :: ites ! max harmonic degree used for Earth field ! name of harmonic coef/const files CHARACTER*(50), PUBLIC :: modfile ! data on the planet we are orbiting around ! GM_{planet}, R_{planet}; REAL(KIND=dkind), PUBLIC :: gm_planet, r_planet ! spherical harmonics coefficients, non normalized REAL(KIND=dkind), PUBLIC :: harmco(narmx) ! public routines PUBLIC :: parm10, geopot, numcoe CONTAINS SUBROUTINE geopot(ites) ! GEOPOT: input of geopotential coefficients ! name of potential model file: modfile (shared data) ! max degree, normalization INTEGER, INTENT(IN) :: ites ! harmonic coefficients, normalized if armno=T ! are placed in the vector harmco ! indexes, scalar temporaries INTEGER :: lll, ind, n ! harmonic indexes and coefficients INTEGER :: l,m,i, ncoef REAL(KIND=dkind) :: c, s ! unit for coefficients file INTEGER :: iun ! name of potential model file CHARACTER (LEN=60) :: modelfile ! dimensions check ncoef=numcoe(ites) IF(ncoef.gt.narmx)THEN write(*,*) 'geopot: ncoef =',ncoef,' narmx=',narmx STOP ENDIF ! name of the geopotential coefficients files modelfile=modfile call rmsp(modelfile,lll) ! assign planet whose harmonic expansion is used ! place radius, mass of planet in spher\_harm.mod ! should be consistent with gravity model r_planet=reau gm_planet=gmearth ! read coefficients in [modfile].coe file CALL filopn(iun,modelfile(1:lll)//'.coe','old') ! default coefficients are zero DO n=1,narmx read(iun,*,end=2)l,m,c, s ! selection of the harmonics up to degree l=ites ! that is the first numcoe(ites) coefficients if(l.gt.ites) CYCLE ind=indcoe(l,m,1) harmco(ind)=c/fnspar(l,m) IF(m.ne.0)THEN ind=indcoe(l,m,0) harmco(ind)=s/fnspar(l,m) ENDIF ENDDO 2 CALL filclo(iun,' ') END SUBROUTINE geopot SUBROUTINE parm10 (x,ites,arm,partials) ! parm: vers. orbit12, revised October 1999 ! computes accelerations due to Mercury spherical ! harmonics up to degree $l=ites > 1$. ! x= body fixed cartesian coordinates ! ites = maximum degree l ! arm(k,1,lmi)= acceleration due to harmonic lmi, k=1,3 ! arm(k,j+1,lmi)= partial derivative of component k of accel. with resp ! component j of x, only if partials=T} ! Warning: to be optimized for performance ! INPUT REAL(KIND=dkind), DIMENSION(3), INTENT (IN) :: x INTEGER, INTENT (IN) :: ites LOGICAL, INTENT(IN) :: partials ! OUTPUT REAL(KIND=dkind), DIMENSION(3,4,narmx), INTENT(OUT) :: arm ! INCLUDE 'planet.h' moved to module shared data REAL(KIND=dkind), DIMENSION(narmx) :: pt,p REAL(KIND=dkind), DIMENSION(itesx) :: s,c,alfa REAL(KIND=dkind), DIMENSION(3) :: dvy REAL(KIND=dkind), DIMENSION(3,3) :: ddv, dy REAL(KIND=dkind), DIMENSION(3,3,3) :: t REAL(KIND=dkind) :: x12,x22,x32,r,r2,r3,rl,rt,rxy,rxy2,rxz2,ryz2, & & den2,rxyq REAL(KIND=dkind) :: d,sf,cf,tgf,cl,sl,pp INTEGER :: ii,lmi ! indexes INTEGER :: i,j,k,l,m,n,nn SAVE ! Warning: dipole terms do not exist if(ites.le.1)return ! Polar coordinates $(r, \phi, \lambda)$, radius, latitude, longitude: ! \[ \sin(\phi)=\frac{z}{r} \quad \cos(\phi)=\frac{x^2+y^2}{r} \] ! \[ \cos(\lambda)=\frac{x}{x^2+y^2} \quad \sin(\lambda)=\frac{y}{x^2+y^2} \] rt=r_planet rxy2=x(1)**2+x(2)**2 r2=rxy2+x(3)**2 r=dsqrt(r2) rxy=dsqrt(rxy2) sf=x(3)/r cf=rxy/r tgf=x(3)/rxy cl=x(1)/rxy sl=x(2)/rxy ! If partials have to be computed, the tensor t(i,j,k), j.ge.k, ! of the second derivatives of $(r, \phi, \lambda)$ with respect ! to $(x,y,z)$ has to be computed if (partials)then r3=r2*r x12=x(1)**2 x22=x(2)**2 x32=x(3)**2 rxz2=r2-x22 ryz2=r2-x12 den2=r2*r2*rxy**3 rxyq=rxy2**2 ! \[ \der2parz{r}{(xyz)}=\frac{1}{r^3}\begin{pmatrix} ! y^2+z^2 & -xy & -xz \cr ! -xy & x^2 + z^2 & -yz \cr ! -xz & -yz & x^2+y^2 \cr ! \end{pmatrix} \] t(1,1,1)=ryz2/r3 t(1,2,1)=-x(1)*x(2)/r3 t(1,3,1)=-x(1)*x(3)/r3 t(1,2,2)=rxz2/r3 t(1,3,2)=-x(2)*x(3)/r3 t(1,3,3)=rxy2/r3 ! \[ \der2parz{\phi}{(xyz)}=\frac{1}{r^4(x^2+y^2)^{3/2}} \times \] ! \[ \times \begin{pmatrix} ! zx^2(2x^2+y^2)-zy^2(z^2+y^2) & ! idem & idem\cr xyz[3x^2+3y^2+z^2] & ! zy^2(2y^2+x^2)-zx^2(z^2+x^2) & idem\cr ! x[z^2-(x^2+y^2)](x^2+y^2) & ! y[z^2-(x^2+y^2)](x^2+y^2) & -2z(x^2+y^2)^2 \cr ! \end{pmatrix} \] t(2,1,1)=x(3)*(x12*(x12+rxy2)-x22*ryz2)/den2 t(2,2,1)=x(1)*x(2)*x(3)*(x12+2.d0*rxy2+ryz2)/den2 t(2,3,1)=x(1)*(x32-rxy2)*rxy2/den2 t(2,2,2)=x(3)*(x22*(x22+rxy2)-x12*rxz2)/den2 t(2,3,2)=x(2)*(x32-rxy2)*rxy2/den2 t(2,3,3)=-2*x(3)*rxyq/den2 ! \[ \der2parz{\lambda}{(xyz)}= \frac{1}{(x^2+y^2)^2}\begin{pmatrix} ! 2xy & y^2-x^2 & 0 \cr y^2-x^2 & -2xy & 0 \cr 0 & 0 & 0 \cr ! \end{pmatrix} \] t(3,1,1)=2*x(1)*x(2)/rxyq t(3,2,1)=(x22-x12)/rxyq t(3,3,1)=0.d0 t(3,2,2)=-2*x(1)*x(2)/rxyq t(3,3,2)=0.d0 t(3,3,3)=0.d0 endif ! $dy(i,j)=$ jacobian matrix of the polar coordinates ! with respect to the cartesian ones ! \[ dy=\derparz{(r,\phi,\lambda)}{(xyz)} \quad ! \derparz{r}{x}=\frac{x}{r} \, \, etc. \] ! \[ \derparz{\phi}{x}=\frac{-xz}{r^2 (x^2+y^2)} \quad ! \derparz{\phi}{y}=\frac{-yz}{r^2 (x^2+y^2)} \quad ! \derparz{\phi}{z}=\frac{x^2+y^2}{r^2} \] ! \[ \derparz{\lambda}{x}=\frac{-y}{x^2+y^2} \quad ! \derparz \lambda y=\frac{x}{x^2+y^2} \] dy(1,1)=x(1)/r dy(1,2)=x(2)/r dy(1,3)=x(3)/r dy(2,1)=-x(1)*x(3)/(r2*rxy) dy(2,2)=-x(2)*x(3)/(r2*rxy) dy(2,3)=rxy/r2 dy(3,1)=-x(2)/rxy2 dy(3,2)=x(1)/rxy2 dy(3,3)=0.d0 ! Legendre associated functions ! \[ p(lmi)=P_{l,m}(\sin\phi) \ \ for\ i=1 \] ! \[ pt(lmi)=P_{l,m}(\sin\phi)\begin{cases} ! \cos(m\lambda) \ \ for\ i=1\cr ! \sin(m\lambda) \ \ for\ i=0,m\neq 0 ! \end{cases} \] ii=(ites+1)**2 call legit(sf,cf,sl,cl,p,pt,s,c,ii,ites) ! \[ \alpha(l)=\frac {GM}{r} \frac{R_\oplus^l}{r^l} \] alfa(1)=gm_planet*rt/r2 ! \[ U=\sum_{l,m,i} J_{l,m,i} U_{l,m,i} \quad ! U_{l,m,i}=\alpha(l)P_{lm}(\sin\phi) ! \begin{cases} \cos(m\lambda)\, \, i=1\cr \sin(m\lambda) \, \, i=0, m\neq 0 ! \end{cases} \] ! loop on degree l do 31 l=2,ites alfa(l)=alfa(l-1)*rt/r rl=-(l+1)/r ! loop on order m do 5 m=0,l ! sine and cosine terms do 51 i=0,1 lmi=l*l+2*m+i if(m.eq.0.and.i.eq.0)then ! for $m=0$ cosine only CYCLE endif ! gradient of the potential in polar coordinates ! \[ \derparz{U_{l,m,i}}{r}= \frac{-(l+1)}{r} \alpha(l)\, P_{lm} ! \begin{cases} \cos(m\lambda)\, \, i=1\cr \sin(m\lambda\, \, i=0 ! \end{cases} \] dvy(1)=rl*alfa(l)*pt(lmi) ! \[ \derparz{U_{l,m,i}}{\lambda}= m \alpha(l)\, P_{lm} ! \begin{cases} -\sin(m\lambda)\, \, i=1 \cr \cos(m\lambda)\, \, i=0 ! \end{cases} \] dvy(3)=(-1)**i*m*alfa(l)*pt(lmi+1-2*i) ! \[ \derparz{U_{l,m,i}}{\phi}= \alpha(l) P'_{lm} \cos\phi ! \begin{cases} \cos(m\lambda)\, \, i=1\cr \sin(m\lambda)\, \, i=0 ! \end{cases} \] ! \[ P'_{lm}\cos\phi=P_{l,m+1}-m \tan\phi\, P_{lm} \] ! Warning: $P_{l,l+1}=0$ if(m.eq.l)then dvy(2)=alfa(l)*(-m*tgf*pt(lmi)) elseif(m.eq.0)then dvy(2)=alfa(l)*p(lmi+2) else pp=p(l*l+2*m+3)*(s(m)*(1-i)+c(m)*i) dvy(2)=alfa(l)*(pp-m*tgf*pt(lmi)) endif ! gradient of the potential in cartesian coordinates ! \[ \derparz{U}{(xyz)}= ! \derparz{U}{(r\phi\lambda)} \derparz{(r\phi\lambda)}{(xyz)} \] DO j=1,3 d=0.d0 DO k=1,3 d=d+dvy(k)*dy(k,j) ENDDO arm(j,1,lmi)=d ENDDO if(partials)then ! ddv(i,j): matrix of second derivatives of U ! \[ \der2parz{U_{l,m,i}}{r} = \frac{(l+1)(l+2)}{r^2}\alpha(l)\, ! P_{lm} \begin{cases} \cos(m\lambda)\cr ! \sin(m\lambda) ! \end{cases} \] ddv(1,1)=-(l+2)*dvy(1)/r ! \[ \frac{\partial^2 U_{l,m,i}}{\partial r \partial\phi}= ! \frac{-(l+1)}r \alpha(l) P'_{lm} \begin{cases} \cos(m\lambda)\cr ! \sin(m\lambda) ! \end{cases}= \frac{-(l+1)}{r}\derparz{U_{l,m,i}}{\phi} \] ddv(1,2)=-(l+1)*dvy(2)/r ddv(2,1)=ddv(1,2) ! \[ \frac{\partial^2 U_{l,m,i}}{\partial r \partial\lambda}= ! \frac{-(l+1)}{r} \alpha(l) m\, P_{lm}\begin{cases} ! -\sin(m\lambda) \cr \cos(m\lambda) ! \end{cases} ! = \frac{-(l+1)}r \derparz {U_{l,m,i}}{\lambda} \] ddv(1,3)=-(l+1)*dvy(3)/r ddv(3,1)=ddv(1,3) ! \[ \der2parz{U_{l,m,i}}{\phi}= \alpha(l) \der2parz{P_{lm}(\sin\phi)}{ ! \begin{cases} \cos(m\lambda) \cr \sin(m\lambda) ! \end{cases}} \] ! \[ \der2parz{P_{lm}(\sin\phi)}{\phi}=\tan\phi P_{l,m+1} + ! [m^2 \sec^2\phi -m\tan^2\phi-l(l+1)] P_{lm}(\sin\phi) \] ddv(2,2)=tgf*dvy(2)+alfa(l)*pt(lmi)*((m/cf)**2-l*(l+1)) ! \[ \frac{\partial^2 U_{l,m,i}}{\partial\lambda \partial\phi}= ! m \alpha(l) P'_{lm} \cos\phi \begin{cases} ! -\sin(m\lambda) \, \, i=1\cr ! \cos(m\lambda) \, \, i=0 ! \end{cases} \] ! \[ P'_{lm}\cos\phi=P_{l,m+1}-m \tan\phi\ P_{lm} \] ! Warning: $P_{l,l+1}=0$ if(m.eq.l)then ddv(2,3)=m*alfa(l)*(-tgf*m*pt(lmi+1-2*i))*(-1)**i elseif(m.eq.0)then ddv(2,3)=0.d0 else pp=p(l*l+2*m+3)*((1-i)*c(m)+i*s(m)) ddv(2,3)=m*alfa(l)*(pp- & & tgf*m*pt(lmi+1-2*i))*(-1)**i endif ddv(3,2)=ddv(2,3) ! \[ \der2parz U\lambda= -m^2 \alpha(l)\ P_{lm} \begin{cases} \cos(m\lambda)\, ! \sin(m\lambda) ! \end{cases} \] ddv(3,3)=-m*m*alfa(l)*pt(lmi) ! matrix of partials of the accelerations with respect to coordinates ! \[ \der2parz{U}{(xyz)}=\bigl[\derparz{(r\phi\lambda)}{(xyz)}\bigl]^T ! \der2parz{U}{(r\phi\lambda)} \derparz{(r\phi\lambda)}{(xyz)} + ! \derparz{U}{(r\phi\lambda)} \der2parz{(r\phi\lambda)}{(xyz)} \] DO k=1,3 DO j=1,k d=0.d0 DO n=1,3 d=d+t(n,k,j)*dvy(n) DO nn=1,3 d=d+ddv(n,nn)*dy(n,k)*dy(nn,j) ENDDO ENDDO arm(j,k+1,lmi)=d arm(k,j+1,lmi)=d ENDDO ENDDO endif ! end loop on i=0,1 51 ENDDO ! end loop on m=0,l 5 ENDDO ! end loop on l=2,ites 31 ENDDO END subroutine parm10 ! \bighorline SUBROUTINE legit(sf,cf,sl,cl,p,pt,s,c,ii,i) !\comments ! {\bf legit} ! recursive computation of Legendre polynomials ! and Legendre associated functions, unnormalised, ! up to degree $l=i$ ! {\obeylines ! INPUT: ! $sf=\sin\phi\quad cf=\cos\phi$ ! $sl=\sin\lambda\quad cf=\cos\lambda $ ! OUTPUT: ! p,pt=vectors of length $ii=(i+1)^2$: ! $p(l^2+2m+1)=P_{l,m}(\sin\phi)$ Legendre associated function ! $pt(l^2+2m+1)=P_{l,m}(\sin\phi)\cos(m\lambda)$ ! $pt(l^2+2m)=P_{l,m}(\sin\phi)\sin(m\lambda) \ \ \ m\neq 0$} ! \interface ! INPUT REAL(KIND=dkind), INTENT (IN) :: sf,sl,cf,cl INTEGER, INTENT (IN) :: i,ii ! OUTPUT REAL(KIND=dkind), DIMENSION(i), INTENT (OUT) :: c,s REAL(KIND=dkind), DIMENSION(ii), INTENT (OUT) :: p,pt ! \endint INTEGER :: j ! indexes INTEGER :: l,m SAVE ! ! initialization of recursion ! \[ P_{00}(\sin\phi)=1 \quad P_{10}(\sin\phi)=\sin\phi \quad ! P_{11}(\sin\phi)=\cos\phi \] p(1)=1.d0 p(2)=sf p(4)=cf ! \[ c(m)=\cos(m\lambda) \quad s(m)=\sin(m\lambda) \] c(1)=cl s(1)=sl ! recursion formulae: DO l=2,i ! for Legendre polynomials ($m=0$): ! \[ P_{l}(\sin\phi)=\frac 1l [(2l-1)\sin\phi P_{l-1}-(l-1)P_{l-2} \] p(l*l+1)=((2*l-1)*sf*p((l-1)*(l-1)+1)-(l-1)*p((l-2)*(l-2)+1))/l ! for $m=l$: ! \[ P_{ll}(\sin\phi)= (2l-1)\cos\phi P_{l-1,l-1} \] p((l+1)*(l+1))=(2*l-1)*cf*p(l*l) ! for $m=l-1$: ! \[ P_{l,l-1}(\sin\phi)=(2l-1)\cos\phi P_{l-1,l-2} \] p((l+1)*(l+1)-2)=(2*l-1)*cf*p(l*l-2) ! for $1\leq m \leq l-2\ ,\ l\geq 2$: ! \[ P_{l,m}(\sin\phi)=P_{l-2,m}+(2l-1)P_{l-1,m-1}\cos\phi \] j=l-2 IF(j.ge.0)THEN DO m=1,j p(l*l+2*m+1)=p((l-2)*(l-2)+2*m+1)+(2*l-1)*cf & & *p((l-1)*(l-1)+2*m-1) ENDDO ENDIF ! addition formulae for sines and cosines of multiples of longitude s(l)=s(l-1)*cl+c(l-1)*sl c(l)=c(l-1)*cl-s(l-1)*sl ENDDO DO l=1,i ! zonal case pt(l*l+1)=p(l*l+1) ! sectorial and tesseral cases DO m=1,l j=l*l+2*m pt(j)=p(j+1)*s(m) pt(j+1)=p(j+1)*c(m) ENDDO ENDDO END SUBROUTINE legit ! \bighorline INTEGER FUNCTION indcoe(l,m,i) ! \comments ! \begin{verbatim} ! Harmonics addressing routines: $lmi=l^2+2m+i$ ! Warning: list of coefficients begins with $\Delta C_{00}$; ! space is left even for the dipole terms, never used ! \end{verbatim} ! {\bf indcoe}: address of harmonics (l,m,i) ! \interface INTEGER, INTENT (IN) :: i,l,m ! \endint if(l.lt.0)stop' **** indcoe: l < 0 ****' if(m.lt.0)stop' **** indcoe: m < 0 ****' if(m.gt.l)stop' **** indcoe: m > l ****' if(m.eq.0)then if(i.ne.1)stop' **** indcoe: i .ne. 1 (m=0) ****' else if(i.ne.1.and.i.ne.0)stop' **** indcoe: i = ? (m>0) ****' endif indcoe=l*l+2*m+i END FUNCTION indcoe ! \bighorline SUBROUTINE namcoe(lmi,l,m,i) ! \comments ! {\bf namcoe} inverse of indcoe ! \interface INTEGER, INTENT (IN) :: lmi INTEGER, INTENT (OUT) :: l,m,i ! \endint REAL(KIND=dkind) :: sl sl=lmi l=dsqrt(sl) ! fix for case of $S_l0$ not existing IF(l*l-lmi.eq.0)l=l-1 m=(lmi-l*l)/2 i=lmi-l*l-2*m END SUBROUTINE namcoe ! \bighorline INTEGER FUNCTION numcoe(ites) ! \comments ! {\bf numcoe} number of coefficients up to degree ites ! \interface INTEGER, INTENT (IN) :: ites ! \endint numcoe=(ites+1)**2 END FUNCTION numcoe ! \bighorline REAL(KIND=dkind) FUNCTION fnspar(l,m) ! \comments ! \begin{verbatim} ! F N S P A R : Factor for the normalization of spherical harmonic ! coefficients ! ! INPUT: ! L : Degree ! M : Order ! OUTPUT: ! FNSPAR: Ratio between normalized and unnormalized coefficient: ! Cnorm(l,m)=fnspar(l,m)*C(l,m) ! \end{verbatim} ! \interface INTEGER, INTENT (IN) :: l,m ! degree, order ! \endint INTEGER :: ne REAL(KIND=dkind) :: f ! indexes INTEGER :: j ! fnspar=1.d0 ne=0 if(m.ne.0)then do j=l-m+1,l+m fnspar=fnspar*j if(fnspar.gt.1.d10)then fnspar=fnspar/1.d10 ne=ne+1 end if enddo fnspar=fnspar/2.d0 end if fnspar=fnspar/(2*l+1) fnspar=sqrt(fnspar) if(ne.ne.0)then f=1.d5**ne fnspar=f*fnspar end if END FUNCTION fnspar END MODULE spher_harm