ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c tass1.7 : new version including hyperion (5 nov 1996) c ************************************************************************* * * from tass1.6 by a. vienne and l. duriez (1995, a&a 297, 588-605) * and 'theory of motion and ephemerides of hyperion', to appear in a&a * e-mail: vienne@gat.univ-lille1.fr * e-mail: duriez@gat.univ-lille1.fr * * positions and velocities of the satellites mimas, enceladus, tethys, * dione, rhea, titan, hyperion and iapetus referred to the center of * saturn and to the mean ecliptic and mean equinox for j2000.0 epoch * ************************************************************************* subroutine posired(dj,is,xyz,vxyz) implicit double precision (a-h,o-z) dimension xyz(3),vxyz(3),elem(6),dlo(8) if(is.eq.7) then call elemhyp(dj,elem) else call calclon(dj,dlo) call calcelem(dj,is,elem,dlo) end if call edered(elem,xyz,vxyz,is) return end subroutine calcelem(dj,is,elem,dlo) implicit double precision (a-h,o-z) dimension elem(6) parameter(ntmx=250) dimension series(ntmx,3,4,8),ntr(5,8),al0(8),an0(8) dimension iks(ntmx,4,8,8) dimension dlo(8) common /taseries/ series,ntr,al0,an0,iks t=(dj-2444240.d0)/365.25d0 s=0 do i=1,ntr(1,is) phas=series(i,2,1,is) do jk=1,8 phas=phas+iks(i,1,is,jk)*dlo(jk) end do s=s+series(i,1,1,is)*cos(phas+t*series(i,3,1,is)) enddo elem(1)=s s=dlo(is)+al0(is) do i=ntr(5,is)+1,ntr(2,is) phas=series(i,2,2,is) do jk=1,8 phas=phas+iks(i,2,is,jk)*dlo(jk) end do s=s+series(i,1,2,is)*sin(phas+t*series(i,3,2,is)) enddo s=s+an0(is)*t cs=cos(s) sn=sin(s) elem(2)=datan2(sn,cs) s1=0 s2=0 do i=1,ntr(3,is) phas=series(i,2,3,is) do jk=1,8 phas=phas+iks(i,3,is,jk)*dlo(jk) end do s1=s1+series(i,1,3,is)*cos(phas+t*series(i,3,3,is)) s2=s2+series(i,1,3,is)*sin(phas+t*series(i,3,3,is)) enddo elem(3)=s1 elem(4)=s2 s1=0 s2=0 do i=1,ntr(4,is) phas=series(i,2,4,is) do jk=1,8 phas=phas+iks(i,4,is,jk)*dlo(jk) end do s1=s1+series(i,1,4,is)*cos(phas+t*series(i,3,4,is)) s2=s2+series(i,1,4,is)*sin(phas+t*series(i,3,4,is)) enddo elem(5)=s1 elem(6)=s2 return end subroutine calclon(dj,dlo) implicit double precision (a-h,o-z) parameter(ntmx=250) dimension series(ntmx,3,4,8),ntr(5,8),al0(8),an0(8) dimension iks(ntmx,4,8,8) dimension dlo(8) common /taseries/ series,ntr,al0,an0,iks t=(dj-2444240.d0)/365.25d0 do is=1,8 if (is.ne.7) then s=0 do i=1,ntr(5,is) s=s+series(i,1,2,is)*sin(series(i,2,2,is)+t*series(i,3,2,is)) end do dlo(is)=s else dlo(is)=0.d0 end if end do return end subroutine lecser(nomf,icrt) implicit double precision (a-h,o-z) parameter(ntmx=250) dimension series(ntmx,3,4,8),ntr(5,8),al0(8),an0(8) dimension iks(ntmx,4,8,8),ik(8) character nomf*(*) dimension am(9),aam(9) dimension tam(9),tmas(9) common /taseries/ series,ntr,al0,an0,iks common /myrd/ aam common /edre/ aia,oma,tmas,gk1 common /trigo/ pi,radsdg save radsdg=datan(1.d0)/45.d0 pi=180.d0*radsdg open(15,file=nomf,status='old') read(15,*) gk read(15,*) tas gk1=(gk*365.25d0)**2.d0/tas read(15,*) aia,oma aia=aia*radsdg oma=oma*radsdg read(15,*) (tam(i),i=1,9) do i=1,9 tmas(i)=1.d0/tam(i) end do read(15,*) (am(i),i=1,9) do i=1,9 aam(i)=am(i)*365.25d0 enddo call dzer(series,ntmx*3*4*8) call nzer(ntr,5*8) call dzer(al0,8) call dzer(an0,8) 10 read(15,*,end=11) is,ieq if(is.eq.7) goto 17 if(ieq.eq.2)then read(15,*) nt, al0(is),an0(is) endif kt=0 isaut=0 9 read(15,*) nt,a1,a2,a3,ik if (nt.lt.9998) then if(kt.eq.ntmx)then write(*,*)'ntmx trop petit' stop endif if (isaut.eq.1) goto 9 kt=kt+1 series(kt,1,ieq,is)=a1 series(kt,2,ieq,is)=a2 series(kt,3,ieq,is)=a3 do js=1,8 iks(kt,ieq,is,js)=ik(js) end do goto 9 else if (nt.eq.9998) then if (icrt.eq.1) isaut=1 if (ieq.eq.2) ntr(5,is)=kt goto 9 end if if(nt.eq.9999) then if ((ieq.eq.2).and.(ntr(5,is).eq.0)) ntr(5,is)=kt ntr(ieq,is)=kt goto 10 end if end if 17 call lithyp(15) 11 close(15) return end subroutine edered(elem,xyz,vxyz,isat) implicit double precision (a-h,o-z) parameter(eps=1.d-14) dimension elem(6),xyz(3),vxyz(3) dimension aam(9),tmas(9) dimension xyz2(3),vxyz2(3) common /myrd/ aam common /edre/ aia,oma,tmas,gk1 amo=aam(isat)*(1.d0+elem(1)) rmu=gk1*(1.d0+tmas(isat)) dga=(rmu/(amo*amo))**(1.d0/3.d0) rl=elem(2) rk=elem(3) rh=elem(4) fle=rl-rk*dsin(rl)+rh*dcos(rl) 20 continue cf=dcos(fle) sf=dsin(fle) corf=(rl-fle+rk*sf-rh*cf)/(1-rk*cf-rh*sf) fle=fle+corf if (dabs(corf).ge.eps) goto 20 cf=dcos(fle) sf=dsin(fle) dlf=-rk*sf+rh*cf rsam1=-rk*cf-rh*sf asr=1.d0/(1.d0+rsam1) phi=dsqrt(1.d0-rk*rk-rh*rh) psi=1.d0/(1.d0+phi) x1=dga*(cf-rk-psi*rh*dlf) y1=dga*(sf-rh+psi*rk*dlf) vx1=amo*asr*dga*(-sf-psi*rh*rsam1) vy1=amo*asr*dga*( cf+psi*rk*rsam1) dwho=2.d0*dsqrt(1.d0-elem(6)*elem(6)-elem(5)*elem(5)) rtp=1.d0-2.d0*elem(6)*elem(6) rtq=1.d0-2.d0*elem(5)*elem(5) rdg=2.d0*elem(6)*elem(5) xyz2(1) = x1 * rtp + y1 * rdg xyz2(2) = x1 * rdg + y1 * rtq xyz2(3) = ( - x1 * elem(6) + y1 * elem(5) ) * dwho vxyz2(1)= vx1 * rtp + vy1 * rdg vxyz2(2)= vx1 * rdg + vy1 * rtq vxyz2(3)= ( -vx1 * elem(6) + vy1 * elem(5) ) * dwho ci=dcos(aia) si=dsin(aia) co=dcos(oma) so=dsin(oma) xyz(1)=co * xyz2(1) - so*ci * xyz2(2) + so*si * xyz2(3) xyz(2)=so * xyz2(1) + co*ci * xyz2(2) - co*si * xyz2(3) xyz(3)= si * xyz2(2) + ci * xyz2(3) vxyz(1)=co * vxyz2(1) - so*ci * vxyz2(2) + so*si * vxyz2(3) vxyz(2)=so * vxyz2(1) + co*ci * vxyz2(2) - co*si * vxyz2(3) vxyz(3)= si * vxyz2(2) + ci * vxyz2(3) return end subroutine nzer(itab,n) integer itab(n) do 10 i=1,n itab(i)=0 10 continue return end subroutine dzer(rtab,n) real*8 rtab(n) do 10 i=1,n rtab(i)=0.d0 10 continue return end subroutine lithyp(nf) implicit real*8(a-h,o-z) parameter(ntp=120, ntq=240, ntz=200, ntzt=65) dimension serp(ntp),serq(ntq),serz(ntz),serzt(ntzt) dimension fap(ntp),faq(ntq),faz(ntz),fazt(ntzt) dimension frp(ntp),frq(ntq),frz(ntz),frzt(ntzt) common /serhyp/nbtp,nbtq,nbtz,nbtzt,t0,cstp,cstq,amm7, & serp,fap,frp, & serq,faq,frq, & serz,faz,frz, & serzt,fazt,frzt save read(nf,'(f15.2)') t0 read(nf,'(d25.16)') amm7 read(nf,'(i5)') nbtp read(nf,'(d25.16)') cstp do i=1,nbtp read(nf,'(3d25.16)') serp(i),fap(i),frp(i) end do read(nf,'(i5)') nbtq read(nf,'(d25.16)') cstq do i=1,nbtq read(nf,'(3d25.16)') serq(i),faq(i),frq(i) end do read(nf,'(i5)') nbtz do i=1,nbtz read(nf,'(3d25.16)') serz(i),faz(i),frz(i) end do read(nf,'(i5)') nbtzt do i=1,nbtzt read(nf,'(3d25.16)') serzt(i),fazt(i),frzt(i) end do return end subroutine elemhyp(dj,elem) implicit real*8(a-h,o-z) dimension elem(6) parameter(pi2=2*3.14159265358979323846d0) parameter(ntp=120, ntq=240, ntz=200, ntzt=65) dimension serp(ntp),serq(ntq),serz(ntz),serzt(ntzt) dimension fap(ntp),faq(ntq),faz(ntz),fazt(ntzt) dimension frp(ntp),frq(ntq),frz(ntz),frzt(ntzt) common /serhyp/nbtp,nbtq,nbtz,nbtzt,t0,cstp,cstq,amm7, & serp,fap,frp, & serq,faq,frq, & serz,faz,frz, & serzt,fazt,frzt t = dj - t0 p=cstp do i=1,nbtp wt=t*frp(i)+fap(i) p = p + serp(i)*dcos(wt) end do q=cstq do i=1,nbtq wt=t*frq(i)+faq(i) q = q + serq(i)*dsin(wt) end do zr=0.d0 zi=0.d0 do i=1,nbtz wt=t*frz(i)+faz(i) zr = zr + serz(i)*dcos(wt) zi = zi + serz(i)*dsin(wt) end do ztr=0.d0 zti=0.d0 do i=1,nbtzt wt=t*frzt(i)+fazt(i) ztr = ztr + serzt(i)*dcos(wt) zti = zti + serzt(i)*dsin(wt) end do vl=mod(amm7*t+q,pi2) if(vl.lt.0.d0) vl=vl+pi2 elem(1)=p elem(2)=vl elem(3)=zr elem(4)=zi elem(5)=ztr elem(6)=zti return end