1subroutine libration(jd,RA,Dec,xl,b) 2 3! Compute optical libration of moon at jd: that is, the sub-observer 4! point (xl,b) in selenographic coordinates. RA and Dec are 5! topocentric values. 6 7 implicit real*8 (a-h,o-z) 8 parameter (RADS=0.0174532925199433d0) 9 parameter (TWOPI=6.28318530717959d0) 10 real*8 jd,j2000,mjd,lambda 11 12 j2000=2451545.0d0 13 RA2000=RA 14 Dec2000=Dec 15 year=2000.0d0 + (jd-j2000)/365.25d0 16 mjd=jd-2400000.d0 17 call sla_PRECES('FK5',year,2000.d0,RA2000,Dec2000) 18 call sla_EQECL(RA2000,Dec2000,mjd,lambda,beta) 19 day=jd - j2000 20 t = day / 36525.d0 21 xi = 1.54242 * RADS 22 ft = 93.2720993 + 483202.0175273 * t - .0034029 * t * t 23 b= ft / 360 24 a = 360 * (b - floor(b)) 25 if (a.lt.0.d0) a = 360 + a; 26 f=a/57.2957795131d0 27 omega=sla_dranrm(2.182439196d0 - t*33.7570446085d0 + t*t*3.6236526d-5) 28 w = lambda - omega 29 y = sin(w) * cos(beta) * cos(xi) - sin(beta) * sin(xi) 30 x = cos(w) * cos(beta) 31 a = sla_dranrm(atan2(y, x)) 32 xl = a - f 33 if(xl.lt.-0.25*TWOPI) xl=xl+TWOPI !Fix 'round the back' angles 34 if(xl.gt.0.25*TWOPI) xl=xl-TWOPI 35 b = asin(-sin(w) * cos(beta) * sin(xi) - sin(beta) * cos(xi)) 36 37 return 38end subroutine libration 39