; ; Earthshine Analysis Program ; ; ; AUTHORS: ; C. Titus Brown(1), Jason Maron(1), Edwin Kolbe(1), Steve Koonin(1), ; Jo Bruls(2), Glenn Eychaner(2), Hal Zirin(2) ; ; (1) Kellogg Radiation Laboratory, Caltech ; (2) Big Bear Solar Observatory, Caltech ; ; ; reference.pro contains the essential almanac reference data for finding ; lunar parameters, as well as other reference functions like phase data. ; ;************************************************************************* ;**************** Package for finding moon parameters ******************** ;************************************************************************* function dms,degrees,minutes,seconds ; ; Convert degrees, minutes and seconds to decimal degrees ; dms = degrees + minutes/60.d0 + seconds/3600.d0 return,dms end ;************************************************************** pro sidtime,JDT,LONGITUDE,LATITUDE,epsilon_mean,epsilon,eqofe,gmst,gast,lmst,last ; ;* sidtime: 94-04-18 ; ; Computes Greenwich and Local mean and apparent sidereal times ; for a given input of Julian date(s), either scalar or array ; ; Input ; jdt : Julian day (fractional) ; longitude : observer's (geodetic) longitude (+East), radials ; latitude : observer's (geodetic) latitude (+North), radials ; Output ; epsilon_mean : mean obliquity of the ecliptic ; epsilon : true obliquity of the ecliptic ; eqofe : equation of equinoxes ; gmst : Greenwich mean sidereal time ; gast : Greenwich apparent sidereal time ; lmst : Local mean sidereal time ; last : Local apparent sidereal time ; ; Time : Astronomical Almanac, 1994, B4-B6 (1984, S15) ; Explanatory Supplement to the Astronomical Almanac, p. 161 ; J2000.0 = JD 2451545.0 = 2000 Jan 1.5 = 2000 Jan 1, 12:00 noon UT1 ; GMST = Greenwich mean sidereal time ; = Greenwich Hour Angle of mean equinox of date ; GAST = Greenwich apparent sidereal time ; = GHA of true equinox of date ; Eqofe = Equation of equinoxes ; = nutation in longitude multiplied by the cosine of ; the obliquity of the ecliptic ; Very accurate values: Methods of astrodynamics, p. 250-260 ; GAST = GMST + Eqofe ; LMST = Local Mean Sidereal Time = GMST + Eastern longitude ; LAST = Local Apparent Sidereal Time = GAST + Eastern longitude ; GHA = Greenwich Hour Angle of some object ; = GAST - RA ; twopi = 2.d0*!dpi xn = jdt-2451545.d0 ; days since J2000.0 tu = xn/36525.d0 ; Julian centuries since J2000.0 GMST = 67310.54841d0+(876600.d0*3600.d0+8640184.812866d0)*tu+0.093104d0*tu^2-6.2d-6*tu^3 GMST = (GMST/3600.d0) mod 24.d0 GMST = (GMST+24.d0*(GMST lt 0.d0))*15.d0*!Dtor ; (radials) ; ; Mean obliquity of the ecliptic of date, i.e. precession only ; Methods of Astrodynamics, p.14 ; Note: the quadratric term is factor 10 too large here. ; epsilon_mean = -46.845d0 * tu1 - 5.9d-3 * tu1^2. + 1.81d-3 * tu1^3. ; epsilon_mean = dms(23.d0,27.d0,8.26d0 + epsilon_mean) * !Dtor ; ; Explanatory Supplement to the Astronomical Almanac, p. 161 ; Astronomical Almanac, 1984 (1994), B18 ; epsilon_mean = -46.8150d0 * tu - 5.9d-4 * tu^2. + 1.813d-3 * tu^3. epsilon_mean = dms(23.d0,26.d0,21.448d0 + epsilon_mean) * !Dtor ; ; True obliquity also needs to take into account short-period (=nutation) ; Order of magnitude of corrections: 10 arcsecs ; Accurate formulae: ; Methods of Astrodynamics, pp. 250-260 ; Explanatory Supplement, pp. 115-120 ; Simple expression, accurate to about 1 arcsec: ; delta_epsilon = 0.0026d0*cos(125.0d0 - 0.05295d0*xn) +$ 0.0002d0*cos(200.9d0 + 1.97129d0*xn) delta_epsilon = delta_epsilon*!Dtor delta_psi = -0.0048d0*sin(125.0d0 - 0.05295d0*xn) -$ 0.0004d0*sin(200.9d0 + 1.97129d0*xn) delta_psi = delta_psi*!Dtor epsilon = epsilon_mean + delta_epsilon Eqofe = delta_psi * cos(epsilon) ; GAST = GMST + Eqofe GAST = GAST mod twopi LMST = (GMST + longitude) mod twopi LMST = LMST + twopi*(LMST lt 0.) LAST = (GAST + longitude) mod twopi LAST = LAST + twopi*(LAST lt 0.) ; return end ;************************************************** pro topocent,LST,OBS_LAT,OBS_ALT,DEC1,RA1,R1,dec2,ra2,r2 ; ; Topocentric coordinates: apply corrections to RA and declination ; Astronomical Almanac, 1994, D3 ; Note that geocentric observer's coordinates should be used here. ; ;* topocent: 94-04-18 ; ; Input ; lst : Local mean or apparent sidereal time ; obs_lat : Observer's geocentric latitude ; obs_alt : Observer's altitude (meters) ; dec1 : Geocentric declination ; ra1 : Geocentric right ascension ; r1 : Geocentric distance in R_eq ; Output ; dec2 : Topocentric declination ; ra2 : Topocentric right ascension ; r2 : Topocentric distance in R_eq ; R_eq = 6378140.e0 x1 = r1*cos(dec1)*cos(ra1) -$ (1.+obs_alt/R_eq)*cos(obs_lat)*cos(LST) y1 = r1*cos(dec1)*sin(ra1) -$ (1.+obs_alt/R_eq)*cos(obs_lat)*sin(LST) z1 = r1*sin(dec1) -$ (1.+obs_alt/R_eq)*sin(obs_lat) ; r2 = sqrt(x1*x1 + y1*y1 + z1*z1) ; distance (Earth radii R_eq) ra2 = atan(y1,x1) ra2 = ra2 + 2.0*!dpi*(ra2 lt 0.) ; mean topocentric RA (radials) dec2 = asin(z1/r2) ; mean topocentric declination (radials) return end ;************************************************** pro alt_azi,latitude,declination,hour_angle,altitude,azimuth ; ; Compute altitude and azimuth from hour angle and declination ; all angles in radials ; ;* alt_azi: 94-04-18 ; ; Input ; latitude : observer's geodetic latitude (+North) ; declination : object's declination ; hour_angle : object's hour angle ; Output ; altitude : object's altitude above horizon ; azimuth : object's azimuth, clockwise from North ; c1 = -cos(declination)*sin(hour_angle) c2 = sin(declination)*cos(latitude) -$ cos(declination)*cos(hour_angle)*sin(latitude) c3 = sin(declination)*sin(latitude) +$ cos(declination)*cos(hour_angle)*cos(latitude) altitude = float(asin(c3)) ; altitude (radials) azimuth = float(atan(c1,c2)) ; azimuth (radials) [-!pi,!pi] azimuth = azimuth+ 2.*!pi*(azimuth lt 0.) ; [0,2.*!pi] azimuth = azimuth*(abs(altitude) ne 90.) ; for zenith and nadir: azimuth=0 ; ;NOTE the value 90 should be pi/2 if we are working in radiants? epb@March04 return end ;************************************************** pro alt_azi_jm,latitude,declination,hour_angle,altitude,azimuth ; ; Compute altitude and azimuth from hour angle and declination ; all angles in radials ; ;* alt_azi: 94-07-26 ; ; Input ; latitude : observer's geodetic latitude (+North) ; declination : object's declination ; hour_angle : object's hour angle ; Output ; altitude : object's altitude above horizon ; azimuth : object's azimuth, clockwise from North ; c1 = -cos(declination)*sin(hour_angle) c2 = sin(declination)*cos(latitude) -$ cos(declination)*cos(hour_angle)*sin(latitude) c3 = sin(declination)*sin(latitude) +$ cos(declination)*cos(hour_angle)*cos(latitude) bbso_x=cos(latitude) bbso_y=0. bbso_z=sin(latitude) moon_x=cos(declination)*cos(hour_angle) moon_y=cos(declination)*sin(hour_angle) moon_z=sin(declination) coaltitude=bbso_x*moon_x+bbso_y*moon_y+bbso_z*moon_z coaltitude=acos(coaltitude) altitude=3.1415927/2.-coaltitude ;print,latitude ;print,declination,hour_angle ;print,bbso_x,bbso_y,bbso_z ;print,moon_x,moon_y,moon_z ;print,coaltitude,altitude ;altitude = float(asin(c3)) ; altitude (radials) azimuth = float(atan(c1,c2)) ; azimuth (radials) [-!pi,!pi] azimuth = azimuth+ 2.*!pi*(azimuth lt 0.) ; [0,2.*!pi] azimuth = azimuth*(abs(altitude) ne 90.) ; for zenith and nadir: azimuth=0 ; return end ;************************************************* pro refract,obs_lat_geoc,obs_alt,altitude,altitude_refrac ; ; Correct altitudes for atmospheric refraction, ; using standard approximations: Astronomical Almanac, 1994, B62 ; Two different approximations apply, depending on whether the altitude ; is larger than 15 degrees or not. ; For altitudes larger than 15 degrees the first approximation (refrac1) ; is accurate to about 0.1 arcmin, but its accuracy deteriorates rapidly ; for smaller altitudes, especially in abnormal meteorological conditions. ; For lower altitudes the second approximation (refrac2) should be used. ; ;* refract: 94-04-18 ; ; Input ; altitude : uncorrected altitude above horizon ; obs_lat_geoc : observer's geocentric latitude ; obs_alt : observer's height above geodetic ; Output ; altitude_refrac : refraction corrected/apparent altitude above horizon ; ; Earth (atmosphere) & moon characterisics ; R_eq, flat : Almanac, 1994, K5-K6 ; xmu : Methods of Astrodynamics, p.19 ; grav_eq, grav : Astrophysical Formulae, p. 496 ; twopi = 2.0*!pi kb = 1.380622e-23 ; Astrophysical Formulae, p. XI (Preface) amu = 1.660531e-27 ; Allen, AQ, p.13 xmu = 28.97 ; molecular mass of atmosphere grav_eq = 9.7803090 ; equatorial gravitational acceleration (m/s^2) sea_pres = 1013 ; ASSUMPTION for atmos pressure (mbar) at sea level atmos_temp = 273.15 ; ASSUMPTION for temperature ; ; Atmospheric scale height for observer's latitude: Astrophys. Form., p. 496 ; It is not quite clear whether geocentric or geodetic latitude should ; be used, but it does not really matter here ; grav = grav_eq + 5.18552e-2*(sin(obs_lat_geoc))^2. -$ 5.70e-5*(sin(2.*obs_lat_geoc))^2. atmos_scale = kb * atmos_temp /(xmu * amu * grav) atmos_pressure = sea_pres * exp(-obs_alt/atmos_scale) ; reflim = 15.*!Dtor refrac1 = 0.00452d0 * atmos_pressure / $ (atmos_temp * tan(altitude)) * !Dtor ; radials refrac2 = atmos_pressure * (0.1594 + 0.0196 * altitude/!Dtor +$ 2.d-5 * (altitude/!Dtor)^2.) / $ (atmos_temp * (1. + 0.505 * altitude/!Dtor +$ 8.45d-2 * (altitude/!Dtor)^2.)) * !Dtor ; radials altitude_refrac = altitude + (altitude ge reflim)*refrac1 +$ (altitude lt reflim)*refrac2 ; refraction corrected altitude return end ;************************************************** pro precess,jdt,dec1,ra1,dec2,ra2 ; ; Precession corrections to convert positions from jdt to J2000.0 ; Note that these should be applied to mean positions, and that ; nutation corrections can be applied afterwards ; Astronomical Almanac, 1994, B18 ; ;* precess: 94-04-18 ; ; Input ; jdt : Starting Julian date(s) ; dec1 : Declination ; ra1 : Right ascension ; Output ; dec2 : Declination ; ra2 : Right ascension ; tu = (jdt-2451545.0d0)/36525. zetaa = (0.6406161d0*tu + 0.0000839d0*tu^2. + 0.0000050d0*tu^3.) * !Dtor za = (0.6406161d0*tu + 0.0003041d0*tu^2. + 0.0000051d0*tu^3.) * !Dtor thetaa = (0.5567530d0*tu - 0.0001185d0*tu^2. - 0.0000116d0*tu^3.) * !Dtor ; c1 = sin(ra1 - za) * cos(dec1) c2 = cos(ra1 - za) * cos(thetaa) * cos(dec1) + sin(thetaa) * sin(dec1) c3 = -cos(ra1 - za) * sin(thetaa) * cos(dec1) + cos(thetaa) * sin(dec1) ; dec2 = asin(c3) ; mean declination for J2000.0 ra2 = atan(c1,c2) - zetaa ra2 = ra2 + 2.d0*!dpi*(ra2 lt 0.) ; mean RA for J2000.0 return end ;************************************************************** pro moon_almanac, JDIM,alphalib_geoc,betalib_geoc,alphalib_topo,betalib_topo,$ phase_geoc,phase_topo,rem_geoc,rem_topo,$ alt_m,alta_m,cres_tilt,pole_m,asun,bsun,alt_m_jm, $ dec_m=dec_m, ra_m=ra_m, dec_s=dec_s, ra_s=ra_s, res=res, $ phi_s = phi_s, phi_m = phi_m, dect_s = dect_s, dect_m = dect_m, $ phit_s = phit_s, phit_m = phit_m, rat_m = rat_m, xn = xn, tu = tu, gmst = gmst, $ epsil = epsilon, gast = gast, Sun_lon = Sun_lon, Sun_anom = Sun_anom, $ bsun0 = bsun0, asun0 = asun0, dect_cen_m = dect_cen_m, rat_cen_m = rat_cen_m, $ phi_cen = phi_cen, R_ms = R_ms, decp = decp, rap = rap, w2_m = w2_m ; ; Returns, for given julian date(s), the solar and lunar coordinates ; and other relevant data ; ;* moon_almanac: 94-07-14 ;*Revised version of sunmoon, with parameter ;*passing instead of common blocks ; ; Input ; jdim (double) : Julian day ; Output ; alphalib_geoc, betalib_geoc : geocentric libration angles ; alphalib_topo, betalib_topo : topocentric libration angles ; phase_geoc : geocentric lunar phase ; phase_topo : topocentric lunar phase ; rem_geoc : geocentric Moon distance (in Earth radii) ; rem_topo : observer-Moon distance(in Earth radii) ; alt_m : theoretical lunar altitudabove horizon ; alta_m : refraction corrected lunar altitude above horizon ; cres_tilt : tilt of center of crescent wrt lunar equator ; positive toward the North ; pole_m : tilt of lunar pole axis wrt observer's ; horizon normal. Zero=vertical, positive ; clockwise ; asun : Selenographic longitude of the sun. ; bsun : Selenographic lattitude of the sun. ; ; ALL OUTPUT ANGLES AND TIMES IN RADIALS !!! ; ;@timecomm ;@earthcom ;@obscomm ;@suncomm ;@mooncomm ; ; Earth & moon characterisics ; R_eq, flat, AU : Almanac, 1994, K5-K6 ; R_m : Almanac, 1994, D46 ; R_eq = 6.378140d6 ; equatorial radius (meters) R_m = 0.2725 ; lunar radius (R_eq) flat = 0.00335281d0 ; flattening factor of Earth AU = 1.4959787066d+11 ; astronomical unit (meters) twopi = 2.d0*!dpi ; ; Observer's coordinates: Astronomical Almanac, 1994, J6 ; Atmospheric scale height for observer's latitude: Astrophys. Form., p. 496 ; Convert to geodetic coordinates: Methods of Astrodynamics, p. 184 ; readcol, 'telnum', telescope_number, format='f', /SILENT ;telescope_number=1.0 if (telescope_number eq 1) then begin ;BBSO observatory obs_lon_geoc = -116.915 * !Dtor ; radials east obs_lat_geoc = 34.2533 * !Dtor ; radials north obs_alt = 2067.0 ; altitude (meters) endif if (telescope_number eq 2) then begin ;Crimea observatory obs_lon_geoc = 34.05 * !Dtor ; radials east obs_lat_geoc = 44.72 * !Dtor ; radials north obs_alt = 586.0 ; altitude (meters) ;When ES3 is at BBSO then---- ;obs_lon_geoc = -116.915 * !Dtor ; radials east ;obs_lat_geoc = 34.2533 * !Dtor ; radials north ;obs_alt = 2067.0 ; altitude (meters) endif if (telescope_number eq 3) then begin ;When ES3 is at BBSO then---- obs_lon_geoc = -116.915 * !Dtor ; radials east obs_lat_geoc = 34.2533 * !Dtor ; radials north obs_alt = 2067.0 ; altitude (meters) ;else ;obs_lon_geoc = 34.05 * !Dtor ; radials east ;obs_lat_geoc = 44.72 * !Dtor ; radials north ;obs_alt = 586.0 ; altitude (meters) endif if (telescope_number eq 4) then begin ;BBSO observatory new ES$ robotic telescope obs_lon_geoc = -116.915 * !Dtor ; radials east obs_lat_geoc = 34.2533 * !Dtor ; radials north obs_alt = 2067.0 ; altitude (meters) endif ;;; Alaska observatory ;; ;;; Alaska observatory ;; ;;;o_lon = -147.75 * !Dtor ; radials east ;;o_lat = 64.83 * !Dtor ; radials north obs_lon_geod = obs_lon_geoc obs_lat_geod = atan(tan(obs_lat_geoc) / (1.d0 - flat)^2.) ;obs_lat_geod = 0.60096159 = 34.432562*!Dtor ; ; Greenwich and local mean and apparent sidereal times and ; mean and true obliquity of the ecliptic ; sidtime,jdim,obs_lon_geod,obs_lat_geod,epsilon_mean,epsilon,eqofe, $ gmst,gast,lmst,last xn = jdim-2451545.d0 ; days since J2000.0 tu = xn/36525.d0 ; Julian centuries since J2000.0 tu1 = tu + 1.d0 ; Julian centuries since 1900 Jan 0.5 ; ; Sun ; Low precision formulae, Astronomical Almanac, 1994, C24 ; accuracy: 0.01 degree (!) for periods between 1950 and 2050 ; Note: the 1984 almanac gives the constant term in Sun_lon as 280.460, but ; the 1994 version lists the updated 280.466 ; Sun_lon = 280.466d0+0.9856474d0*xn ; mean longitude of Sun Sun_anom = 357.528d0+0.9856003d0*xn ; mean anomaly of Sun Sun_lon = Sun_lon mod 360.d0 ; reduce to interval -360 to 360 Sun_anom = Sun_anom mod 360.d0 Sun_lon = Sun_lon+360.d0*(Sun_lon lt 0.d0) ; reduce to interval 0 to 360 Sun_anom = (Sun_anom +360.d0*(Sun_anom lt 0.d0)) * !Dtor ; radials Sun_lon = Sun_lon+1.915d0*sin(Sun_anom)+2.0d-2*sin(2*Sun_anom) ; ecliptic longitude (degrees) Sun_lon = Sun_lon * !Dtor ; radials xbeta = 0.d0 ; ecliptic latitude (definition) res = 1.00014d0 - 0.01671d0*cos(Sun_anom) -$ 0.00014d0*cos(2.*Sun_anom) ; E-S distance in A.U. sd_s = 0.267/res * !Dtor ; semi-diameter (radials) res = res * AU / R_eq ; E-S distance (R_eq) ; decm_s = asin(sin(Sun_lon)*sin(epsilon_mean)) ; mean declination of date ram_s = atan(cos(epsilon_mean)*sin(Sun_lon),cos(Sun_lon)) ; mean RA of date ram_s = ram_s + twopi*(ram_s lt 0.d0) phim_s = GMST - ram_s ; mean Greenwich Hour Angle phim_s = phim_s + twopi*(phim_s lt 0.d0) dec_s = asin(sin(Sun_lon)*sin(epsilon)) ; true declination (radials) ra_s = atan(cos(epsilon)*sin(Sun_lon),cos(Sun_lon)) ; RA (radials) ra_s = ra_s + twopi*(ra_s lt 0.d0) ; interval [0,2*!pi] phi_s = GAST - ra_s ; Greenwich Hour Angle phi_s = phi_s + twopi*(phi_s lt 0.d0) ; ; Moon ; Astronomical Almanac, 1984, S31 : Recommended values for the direction ; of the North pole of Rotation and the Prime Meridian ; for standard epoch 2000 Jan 1.5 TDB, i.e. JD 2451545.0 ; No need to do range checks! ; E1 = (125.045 - 0.052992*xn)*!Dtor E2 = (249.390 - 0.105984*xn)*!Dtor E3 = (196.694 - 13.012000*xn)*!Dtor E4 = (176.630 + 13.340716*xn)*!Dtor E5 = (358.219 - 0.985600*xn)*!Dtor pra2_m = (270.000 - 3.878*sin(E1) - 0.120*sin(E2) +$ 0.070*sin(E3) - 0.017*sin(E4))*!Dtor pdec2_m = (66.534 + 1.543*cos(E1) + 0.024*cos(E2) -$ 0.028*cos(E3) + 0.007*cos(E4))*!Dtor w2_m = (38.314 + 13.1763581*xn + 3.558*sin(E1) + 0.121*sin(E2) -$ 0.064*sin(E3) + 0.016*sin(E4) + 0.025*sin(E5))*!Dtor ; ; Methods of Astrodynamics, p.11-12 : Position ; Moon = dms(270.d0,26.d0,2.99d0) + $ dms(1336.d0*360.d0+307.d0,52.d0,59.31d0) * tu1 - $ dms(0.d0,0.d0,4.08d0) * tu1^2. + $ dms(0.d0,0.d0,0.0068d0) * tu1^3. Gamma = dms(281.d0,13.d0,15.00d0) + $ dms(0.d0,0.d0,6189.03d0) * tu1 + $ dms(0.d0,0.d0,1.63d0) * tu1^2. + $ dms(0.d0,0.d0,0.012d0) * tu1^3. Gammap = dms(334.d0,19.d0,46.40d0) + $ dms(11.d0*360d0+109.d0,2.d0,2.52d0) * tu1 - $ dms(0.d0,0.d0,37.17d0) * tu1^2. - $ dms(0.d0,0.d0,0.045d0) * tu1^3. Omega = dms(259.d0,10.d0,59.79d0) - $ dms(5.d0*360.d0+134.d0,8.d0,31.23d0) * tu1 + $ dms(0.d0,0.d0,7.48) * tu1^2. + $ dms(0.d0,0.d0,0.008d0) * tu1^3. D = dms(350.d0,44.d0,14.95d0) + $ dms(1236.d0*360.d0 + 307.d0,6.d0,51.18d0) * tu1 - $ dms(0.d0,0.d0,5.17d0) * tu1^2. + $ dms(0.d0,0.d0,0.0068d0) * tu1^3. Moon = Moon * !Dtor ; Mean longitude of the moon, measured in the ? ecliptic from the mean equinox of date to the ; mean ascending node of the lunar orbit and then ; along the orbit Gamma = Gamma * !Dtor ; Sun's mean longitude of perigee Gammap = Gammap * !Dtor ; Mean longitude of the lunar perigee, measured ; as Moon Omega = Omega * !Dtor ; Longitude of the mean ascending node of the ; lunar orbit on the ecliptic, measured ; from the mean equinox of date D = D * !Dtor ; Mean elongation of the Moon from the Sun llc = Moon - Gammap luc = Moon - D lp = luc - Gamma F = Moon - Omega ; ; Moon's longitude ; Moon_lon = 22639.580d0 * sin(llc) -$ 4586.438d0 * sin(llc-2.*D) +$ 2369.899d0 * sin(2.*D) +$ 769.021d0 * sin(2.*llc) -$ 668.944d0 * sin(lp) -$ 411.614d0 * sin(2.*F) -$ 211.658d0 * sin(2.*llc - 2.*D) -$ 206.219 * sin(llc + lp - 2.*D) +$ 191.954d0 * sin(llc + 2.*D) -$ 165.351 * sin(lp - 2.*D) +$ 147.878d0 * sin(llc - lp) -$ 124.785d0 * sin(D) -$ 109.804d0 * sin(llc+lp) -$ 55.174d0 * sin(2.*F - 2.*D) -$ 45.100d0 * sin(llc + 2.*F) +$ 39.532d0 * sin(llc - 2.*F) -$ 38.428d0 * sin(llc - 4.*D) +$ 36.124d0 * sin(3.*llc) -$ 30.773d0 * sin(2.*llc - 4.*D) -$ 28.511d0 * sin(llc - lp - 2.*D) -$ 24.451d0 * sin(lp + 2.*D) Moon_lon = (dms(0.d0,0.d0,Moon_lon)*!Dtor + Moon) mod twopi Moon_lon = Moon_lon + twopi*(Moon_lon lt 0.d0) Moon_lat = 18461.480d0 * sin(F) +$ 1010.180d0 * sin(llc + F) -$ 999.695d0 * sin(F - llc) -$ 623.658d0 * sin(F - 2.*D) +$ 199.485d0 * sin(F + 2.*D - llc) -$ 166.577d0 * sin(llc + F - 2.*D) +$ 117.262d0 * sin(F + 2.*D) +$ 61.913d0 * sin(2.*llc + F) -$ 33.359d0 * sin(F - 2.*D - llc) -$ 31.763d0 * sin(F - 2.*D) -$ 29.689d0 * sin(lp + F - 2.*D) Moon_lat = dms(0.d0,0.d0,Moon_lat)*!Dtor Moon_pi = 3422.5400d0 + 186.5398d0*cos(llc) +$ 34.3117d0 * cos(llc-2.*D) +$ 28.2333d0 * cos(2.*D) +$ 10.1657d0 * cos(2.*llc) +$ 3.0861d0 * cos(llc+2.*D) +$ 1.9202d0 * cos(lp-2.*D) +$ 1.4455d0 * cos(llc+lp-2.*D) +$ 1.1542d0 * cos(llc-lp) -$ 0.9752d0 * cos(D) -$ 0.9502d0 * cos(llc+lp) -$ 0.7136d0 * cos(llc-2.*F) +$ 0.6215d0 * cos(3.*llc) +$ 0.6008d0 * cos(llc-4.*D) Moon_pi = dms(0.d0,0.d0,Moon_pi)*!Dtor ; ; Convert ecliptic longitude and latitude to RA and declination ; for the equinox of date ; c1 = cos(Moon_lat)*cos(Moon_lon) c2 = cos(epsilon)*cos(Moon_lat)*sin(Moon_lon) -$ sin(epsilon)*sin(Moon_lat) c3 = sin(epsilon)*cos(Moon_lat)*sin(Moon_lon) +$ cos(epsilon)*sin(Moon_lat) dec_m = asin(c3) ; true declination of date (radials) ra_m = atan(c2,c1) ra_m = ra_m + twopi*(ra_m lt 0.d0) phi_m = GAST - ra_m phi_m = phi_m + twopi*(phi_m lt 0.d0) ; Greenwich hour angle (radials) ;if phi_m gt !pi then phi_m = phi_m - fix(phi_m/!pi)*twopi ; c1 = cos(Moon_lat)*cos(Moon_lon) c2 = cos(epsilon_mean)*cos(Moon_lat)*sin(Moon_lon) -$ sin(epsilon_mean)*sin(Moon_lat) c3 = sin(epsilon_mean)*cos(Moon_lat)*sin(Moon_lon) +$ cos(epsilon_mean)*sin(Moon_lat) decm_m = asin(c3) ; mean declination of date (radials) ram_m = atan(c2,c1) ; mean RA of date (radials) ram_m = ram_m + twopi*(ram_m lt 0.d0) ; mean, i.e. no nutation phim_m = GMST - ram_m phim_m = phim_m + twopi*(phim_m lt 0.d0) ; Greenwich hour angle (radials) ; rem_geoc = 1./tan(Moon_pi) ; E-M distance (Earth radii) sd_m = 0.2725 * Moon_pi ; semi-diameter of Moon rem = rem_geoc ; ; Topocentric coordinates: apply corrections to RA and declination ; Compute both the actual and the mean topocentric coordinates. ; You need the actual coordinates for phase as well as altitude/azimuth ; computation, but you need the mean coordinates to transform to J2000.0 ; ; Mean topocentric RA and declination for mean equinox of date ; topocent,LMST,obs_lat_geoc,obs_alt,decm_m,ram_m,rem_geoc,dectm_m,ratm_m,rem_topo phitm_m = LMST - ratm_m phitm_m = phitm_m + twopi*(phitm_m lt 0.) ; mean topocentric HA (radials) ; topocent,LMST,obs_lat_geoc,obs_alt,decm_s,ram_s,res,dectm_s,ratm_s,res phitm_s = LMST - ratm_s phitm_s = phitm_s + twopi*(phitm_s lt 0.) ; mean topocentric HA (radials) ; topocent,LAST,obs_lat_geoc,obs_alt,dec_m,ra_m,rem_geoc,dect_m,rat_m,rem_topo phit_m = LAST - rat_m ;print,'last, rat_m, phit_m ' ;print,last,rat_m,phit_m phit_m = phit_m + twopi*(phit_m lt 0.) ; true topocentric HA (radials) ;print,twopi ; topocent,LAST,obs_lat_geoc,obs_alt,dec_s,ra_s,res,dect_s,rat_s,res phit_s = LAST - rat_s phit_s = phit_s + twopi*(phit_s lt 0.) ; true topocentric HA (radials) ; new stuff by QJ in 00sep; here want to get the central point on ; the earth's surface between Sun vector and Moon vector, and ; derive the topocentric Moon vector from this point for ; future use. sun_vect = [cos(dec_s)*cos(ra_s), cos(dec_s)*sin(ra_s), sin(dec_s)] moon_vect = [cos(dec_m)*cos(ra_m), cos(dec_m)*sin(ra_m), sin(dec_m)] cen_vect = sun_vect+moon_vect cen_vect = cen_vect/sqrt(total(cen_vect^2)) dec_cen = asin(cen_vect(2)) ra_cen = atan(cen_vect(1), cen_vect(0)) ra_cen = ra_cen + twopi*(ra_cen lt 0.d0) phi_cen = GAST - ra_cen phi_cen = phi_cen + twopi*(phi_cen lt 0.d0) cen_lat_geoc = dec_cen cen_alt = 0.d0 ; altitude (meters) cen_lat_geod = atan(tan(cen_lat_geoc) / (1.d0 - flat)^2.) topocent, ra_cen, cen_lat_geoc,cen_alt,dec_m,ra_m,rem, dect_cen_m, rat_cen_m ; ; Compute altitude and azimuth for Sun and Moon, ; and also indicate if above the horizon ; Note, that now geodetic coordinates are the appropriate ones. ; Astrophysical Formulae, p. 504 ; alt_azi,obs_lat_geod,dect_m,phit_m,alt_m,azi_m alt_azi,obs_lat_geod,dect_s,phit_s,alt_s,azi_s alt_azi_jm,obs_lat_geod,dect_m,phit_m,alt_m_jm,azi_m ; Compute the geocentric and topocentric lunar phases (phase_geoc and ; phase_topo) from the angle between the geocentric and topocentric ; positions of Sun and Moon. ; To get correct positions for where Sun and Moon are at the zenith, ; it is important that you use geocentric coordinates for the ; equinox of date. ; The actual apparent size/phase of the crescent of the moon ; depends on the position on Earth, so that you should ; use topocentric coordinates to compute that! ; ; Note that the angle defined by acos(S.M) is not the lunar phase, but ; (!pi - phase), therefore: cos(phase) = -dot_sm ; cross_sm = -(SxM)_z, instead of +(SxM)_z, simply because hour angles ; are used instead of RA, which are measured in opposite direction, so that the ; y-components of the S and M vectors get the wrong sign (the coordinate system ; defined by HA and declination is not right-handed). ; One can easily see that phase<0 in case (-cross_sm) gt 0, and vice versa ; dot_sm = cos(dect_s)*cos(phit_s)*cos(dect_m)*cos(phit_m) + $ cos(dect_s)*sin(phit_s)*cos(dect_m)*sin(phit_m) + $ sin(dect_s)*sin(dect_m) ; dot product of Sun and Moon vectors (S.M) phase_topo = acos(-dot_sm) ; topocentric lunar phase (radials) cross_sm = cos(dect_s)*cos(phit_s)*cos(dect_m)*sin(phit_m) - $ cos(dect_s)*sin(phit_s)*cos(dect_m)*cos(phit_m) ; -(SxM)_z phase_topo = (2.*(cross_sm ge 0.)-1.) * phase_topo ; topocentric phase dot_sm = cos(dec_s)*cos(phi_s)*cos(dec_m)*cos(phi_m) + $ cos(dec_s)*sin(phi_s)*cos(dec_m)*sin(phi_m) + $ sin(dec_s)*sin(dec_m) ; dot product of Sun and Moon vectors (S.M) phase_geoc = acos(-dot_sm) ; geocentric lunar phase (radials) cross_sm = cos(dec_s)*cos(phi_s)*cos(dec_m)*sin(phi_m) - $ cos(dec_s)*sin(phi_s)*cos(dec_m)*cos(phi_m) ; -(SxM)_z phase_geoc = (2.*(cross_sm ge 0.)-1.) * phase_geoc ; geocentric phase ; ; Note: instead of writing the full cross product formula, you could have ; used function crossp built into IDl, but for one coordinate this direct ; evaluation is simpler ; refract,obs_lat_geoc,obs_alt,alt_m,alta_m refract,obs_lat_geoc,obs_alt,alt_s,alta_s vis_m = (alta_m gt 0.) ; (in)visible = (0)1 vis_s = (alta_s gt 0.) ; (in)visible = (0)1 ; ; Precess the positions of Sun and Moon to the coordinate system for J2000.0 ; required for plotting it on the stellar map ; precess,jdim,dectm_m,ratm_m,dec2_m,ra2_m precess,jdim,dectm_s,ratm_s,dec2_s,ra2_s ; ; Correct dec2_m, ra2_m and dec2_s, ra2_s for nutation, ; when accurate formulae are available, as indicated above, ; from Methods of Astrodynamics, p. 250-260 ; Note, that for plotting this is not really necessary! ; ; Geocentric libration angles alphalib_geoc, betalib_geoc ; These are the same as the selenographic coordinates of the observer ; Results compare well with old method (seleno) and Astronomical Almanac tables ; decm = dec_m ; carefull: all these may be vectors ram = ra_m decp = pdec2_m rap = pra2_m w = w2_m ;moonvec = [cos(decm)*cos(ram),cos(decm)*sin(ram),sin(decm)] ;earthvec = -moonvec ;polevec = [cos(decp)*cos(rap),cos(decp)*sin(rap),sin(decp)] ;e1p = [-sin(rap),cos(rap),0.*rap] ;e2p = [-sin(decp)*cos(rap),-sin(decp)*sin(rap),cos(decp)] ;e3p = polevec ; E1pp = -cos(decm)*cos(ram)*(-cos(w)*sin(rap) -sin(w)*sin(decp)*cos(rap)) $ -cos(decm)*sin(ram)*( cos(w)*cos(rap) -sin(w)*sin(decp)*sin(rap)) $ -sin(decm)*sin(w)*cos(decp) E2pp = -cos(decm)*cos(ram)*( sin(w)*sin(rap) -cos(w)*sin(decp)*cos(rap)) $ -cos(decm)*sin(ram)*(-sin(w)*cos(rap) -cos(w)*sin(decp)*sin(rap)) $ -sin(decm)*cos(w)*cos(decp) E3pp = -cos(decm)*cos(ram)*cos(decp)*cos(rap) $ -cos(decm)*sin(ram)*cos(decp)*sin(rap) $ -sin(decm)*sin(decp) betalib_geoc = asin(E3pp) alphalib_geoc = atan(E2pp,E1pp) ; ; Topocentric libration angles alphalib_topo, betalib_topo ; These are the same as the selenographic coordinates of the observer ; Results compare well with old method (seleno) and Astronomical Almanac tables ; decm = dect_m ; carefull: all these may be vectors ram = rat_m decp = pdec2_m rap = pra2_m w = w2_m E1pp = -cos(decm)*cos(ram)*(-cos(w)*sin(rap) -sin(w)*sin(decp)*cos(rap)) $ -cos(decm)*sin(ram)*( cos(w)*cos(rap) -sin(w)*sin(decp)*sin(rap)) $ -sin(decm)*sin(w)*cos(decp) E2pp = -cos(decm)*cos(ram)*( sin(w)*sin(rap) -cos(w)*sin(decp)*cos(rap)) $ -cos(decm)*sin(ram)*(-sin(w)*cos(rap) -cos(w)*sin(decp)*sin(rap)) $ -sin(decm)*cos(w)*cos(decp) E3pp = -cos(decm)*cos(ram)*cos(decp)*cos(rap) $ -cos(decm)*sin(ram)*cos(decp)*sin(rap) $ -sin(decm)*sin(decp) betalib_topo = asin(E3pp) alphalib_topo = atan(E2pp,E1pp) alib = alphalib_topo blib = betalib_topo ; ; Selenographic coordinates of the Sun ; i.e. direction towards the Sun as seen from the CENTER of the moon xs = res*cos(dect_s)*cos(rat_s) - rem_topo*cos(dect_m)*cos(rat_m) ys = res*cos(dect_s)*sin(rat_s) - rem_topo*cos(dect_m)*sin(rat_m) zs = res*sin(dect_s) - rem_topo*sin(dect_m) R_ms = sqrt(xs*xs + ys*ys +zs*zs) xs = xs/R_ms ys = ys/R_ms zs = zs/R_ms S1pp = xs*(-cos(w)*sin(rap) -sin(w)*sin(decp)*cos(rap)) $ +ys*( cos(w)*cos(rap) -sin(w)*sin(decp)*sin(rap)) $ +zs*sin(w)*cos(decp) S2pp = xs*( sin(w)*sin(rap) -cos(w)*sin(decp)*cos(rap)) $ +ys*(-sin(w)*cos(rap) -cos(w)*sin(decp)*sin(rap)) $ +zs*cos(w)*cos(decp) S3pp = xs*cos(decp)*cos(rap) $ +ys*cos(decp)*sin(rap) $ +zs*sin(decp) bsun = asin(S3pp) asun = atan(S2pp,S1pp) ; these are direction of Sun at zero libration (rap, decp = 0); @QJ ;S1pp0 = ys0 ;S2pp0 = zs0 ;S3pp0 = xs0 ;bsun0 = asin(S3pp0) ;asun0 = atan(S2pp0,S1pp0) ; ; Compute the angle between the lunar pole axis and the observer's ; zenith direction (local 'up'), positive clockwise ; Start with establishing the direction (RA,dec) of local zenith from ; the observer's latitude and local sidereal time (use mean, since you ; are mixing equinoxes anyway ; raz = lmst decz = raz*0. + obs_lat_geod ; make this a vector if raz is a vector ;decp = pdec2_m ; already defined above ;rap = pra2_m ;decm = dect_m ;ram = rat_m ; moonvec = [[cos(decm)*cos(ram)],[cos(decm)*sin(ram)],[sin(decm)]] polevec = [[cos(decp)*cos(rap)],[cos(decp)*sin(rap)],[sin(decp)]] zenithvec = [[cos(decz)*cos(raz)],[cos(decz)*sin(raz)],[sin(decz)]] ; ; Find the projections of polevec and zenithvec onto the plane ; normal to moonvec, by subtracting their projection onto moonvec ; poleproject = total(polevec*moonvec) zenithproject = total(zenithvec*moonvec) polenorm = polevec zenithnorm = zenithvec for i=0,2 do begin polenorm(0,i) = moonvec(*,i)*poleproject zenithnorm(0,i) = moonvec(*,i)*zenithproject endfor polep = polevec - polenorm zenithp = zenithvec - zenithnorm ; ; Renormalize polep and zenithp in order ; to determine the angle between them ; polepnorm = sqrt(total(polep^2.)) zenithpnorm = sqrt(total(zenithp^2.)) polezenith = acos(total(polep*zenithp)/polepnorm/zenithpnorm) polecrosszenith = polevec polecrosszenith(0,0) = polep(*,1)*zenithp(*,2) - polep(*,2)*zenithp(*,1) polecrosszenith(0,1) = polep(*,2)*zenithp(*,0) - polep(*,0)*zenithp(*,2) polecrosszenith(0,2) = polep(*,0)*zenithp(*,1) - polep(*,1)*zenithp(*,0) inprod = total(polecrosszenith*moonvec) pole_m = -polezenith * (2.*(inprod gt 0.) -1.) if (n_elements(pole_m) eq 1) then pole_m = pole_m(0) ; ; Crescent tilt angle wrt to lunar equator. ; Tilt toward the north is positive. ; Start by projecting both the lunar pole axis vector (0,90) and ; the direction (asun,bsun) toward the Sun onto the plane perpendicular ; to the direction toward the observer (alphalib_topo, betalib_topo) ; This is easiest in Cartesian coordinates derived from these selenographic ones. ; ; Construct the vectors and allow for array form of asun and bsun etc. ; polev = [[0*cos(bsun)],[0*cos(bsun)],[0*cos(bsun)+1.]] sunv = [[cos(bsun)*cos(asun)],[cos(bsun)*sin(asun)],[sin(bsun)]] obsv = [[cos(betalib_topo)*cos(alphalib_topo)],$ [cos(betalib_topo)*sin(alphalib_topo)],[sin(betalib_topo)]] ; ; Compute inproducts with observer's direction and subsequently the ; projection onto that direction. ; Subtract the projection vectors to get the vectors in the plane normal ; to the observer's direction ; polevcos = total((polev*obsv)) sunvcos = total((sunv*obsv)) polenorm = obsv sunnorm = obsv for i=0,2 do begin polenorm(0,i) = polev(*,i)*polevcos sunnorm(0,i) = sunv(*,i)*sunvcos endfor polevproj = polev - polenorm sunvproj = sunv - sunnorm ; ; Compute the norm of the projections and their angle ; After that, redefine the angle angle between Sun and lunar north, to ; get the wanted crescent tilt angle with respect to the lunar equator. ; Tilt toward the north being positive angle! ; polepnorm = sqrt(total(polevproj^2.)) sunpnorm = sqrt(total(sunvproj^2.)) cres_tilt = acos(total(polevproj*sunvproj)/polepnorm/sunpnorm) cres_tilt = 0.5*!pi - cres_tilt if (n_elements(cres_tilt) eq 1) then cres_tilt = cres_tilt(0) ; return end ;************************************************************************** ;********** End of lunar parameter package ******************************** ;************************************************************************** ;***************************************************************************** ; Phase_Func_Danjon determines the phase funtion for Danjon's bright-side ; lunar spot by interpolation from his "tableau III" on p.150, ; M.A. Dajon, NOUVELLES RECHERCHES sur la PHOTOMETRIE DE LA LUMIERE ; ET L'ALBEDO DE LA TERRE ; function Phase_Func_Danjon,phase Tableau_in_mag=[0.0 ,0.03,0.06,0.10,0.13,0.17,0.20,0.24,0.27,0.30, $ 0.33,0.36,0.39,0.42,0.45,0.48,0.51,0.53,0.55,0.58, $ 0.60,0.63,0.65,0.68,0.71,0.73,0.75,0.77,0.79,0.81 ] Tableau_in_mag=[Tableau_in_mag, $ 0.83,0.85,0.87,0.89,0.91,0.92,0.93,0.95,0.96,0.98, $ 0.99,1.00,1.01,1.03,1.04,1.05,1.06,1.08,1.09,1.10, $ 1.11,1.12,1.13,1.14,1.15,1.16,1.17,1.18,1.19,1.20 ] Tableau_in_mag=[Tableau_in_mag, $ 1.21,1.22,1.23,1.24,1.25,1.25,1.26,1.27,1.28,1.29, $ 1.29,1.30,1.31,1.32,1.33,1.33,1.34,1.35,1.36,1.37, $ 1.38,1.39,1.40,1.41,1.42,1.43,1.44,1.45,1.46,1.47 ] Tableau_in_mag=[Tableau_in_mag, $ 1.48,1.49,1.50,1.51,1.52,1.53,1.54,1.55,1.56,1.57, $ 1.58,1.59,1.61,1.62,1.63,1.64,1.66,1.67,1.68,1.70, $ 1.71,1.73,1.74,1.76,1.77,1.79,1.80,1.82,1.84,1.86 ] Tableau_in_mag=[Tableau_in_mag, $ 1.88,1.90,1.92,1.94,1.96,1.98,2.01,2.03,2.06,2.08, $ 2.11,2.13,2.16,2.19,2.21,2.24,2.27,2.30,2.33,2.37, $ 2.40,2.43,2.47,2.50,2.54,2.58,2.61,2.66,2.70,2.74 ] Tableau_in_mag=[Tableau_in_mag, $ 2.79,2.84,2.89,2.94,3.00,3.05,3.12,3.18,3.25,3.33, $ 3.41,3.50,3.59,3.69,3.81,3.93] phasedeg=phase*180./3.141593 mag120 = Tableau_in_mag(120) int120 = 10.0^(mag120/(-2.5)) Tab_index=fix(abs(phasedeg)) if (Tab_index le 164) then begin ; interpolation only for phases contained in the Tableau mag=Tableau_in_mag(Tab_index) + (abs(phasedeg)-float(Tab_index)) $ *(Tableau_in_mag(Tab_index+1)-Tableau_in_mag(Tab_index)) int = 10.0^(mag/(-2.5)) ; PFD = int/int120 ; if we want to rescale to 120 PFD = int endif else begin ; otherwise set PFD to a crazy value PFD =-1.0 endelse return,PFD end ; function Phase_Func_Danjon ;***************************************************************************** ; Phase_Func_Opposition_Danjon determines the phase funtion for Danjon's ; bright-side lunar spot by interpolation from his "tableau III" on p.150, ; M.A. Dajon, NOUVELLES RECHERCHES sur la PHOTOMETRIE DE LA LUMIERE ; ET L'ALBEDO DE LA TERRE ; We then correct for the opposition effect by rescaling 20% downwards for ; anything more than 5 deg. away from 0 deg. CTB ; function Phase_Func_Opposition_Danjon,phase pre_opposition = Phase_Func_Danjon(phase) if ((abs(phase)*180.0/3.141592) gt 5.0) then begin val = pre_opposition * .8 endif else begin val = pre_opposition endelse return,val end ; function Phase_Func_Opposition_Danjon