Earthshine blog

Here is the FORTRAN code that predicts the position of the Sun glint given the position of the Moon and the time.

PROGRAM SUNGLINT_MOON
c—————————————————————————
c code that calculates longitude and latitude of the sunglint
c as seen from the Moon
c—————————————————————————
c compile with: gfortran -O3 -fbounds-check sunglintFORTRAN_moon.f -o sung.exe
c or like this: f77 sunglintFORTRAN_moon.f -o sung.exe
c—————————————————————————
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IMPLICIT logical (A-Z)
DOUBLE PRECISION GMT,SANG,DIFF,GLON,GLAT,GLAT1,AERR,GLON1,SLON1
DOUBLE PRECISION SLAT1,DEC1,DLAT,DLON,OLDGLAT1,DEC,RA,ASAT,ZSAT
DOUBLE PRECISION OLDGLON1,GLON2,U1,U2,U3,U4,V1,V2,V3,HANG1,HANG
DOUBLE PRECISION F0,G0,DUDT1,DUDT2,DUDT3,DUDT4,DUDP1,DUDP2,DUDP3
DOUBLE PRECISION DUDP4,HOURANG,ATND,DVDT1,DVDT2,DVDT3,XX1,XX2
DOUBLE PRECISION DVDP1,DVDP2,DVDP3,T1,T2,T3,T4,TEMP1,TEMP2,TEMP3
DOUBLE PRECISION AK,AH,TEMP4,TEMP5,TEMP6,DFDT,DFDP,DGDT,DGDP,JD
DOUBLE PRECISION ALATD,ALOND,ECLAT,ECLON,PI,SLAT,SLON,ALT,AZI
PARAMETER (PI=.3141592654D1)
CHARACTER DATE*15,TIME*8,INSTR*8
INTEGER HOUR,MIN,SEC
c READ IN DATE,TIME(GMT) AND SUB-LUNAR POINT (SLAT,SLON)
PRINT ‘(/5X,A,$)’,’ENTER DATE (DD-MM-YYYY):’
READ ‘(A)’,DATE
PRINT ‘(/5X,A,$)’,’ENTER TIME(GMT) (HH:MM:SS):’
READ ‘(A)’,TIME
PRINT ‘(/5X,A,$)’,’ENTER SUB-LUNAR POSITION (SLAT,SLON)’
READ *,SLAT,SLON
c CONVERT TIME (HH:MM:SS) TO DECIMAL HOURS
WRITE (INSTR,'(A)’) TIME
READ (INSTR,100) HOUR,MIN,SEC
100 FORMAT (3(I2,1X))
GMT=FLOAT(HOUR)+FLOAT(MIN)/.6D2+FLOAT(SEC)/.36D4
c COMPUTE POSITION OF SPECULAR POINT (GLAT,GLON)
CALL ECCS (DATE,GMT,ECLAT,ECLON,JD)
CALL EQCS (DATE,GMT,ECLAT,ECLON,DEC,RA)
CALL NEWRAPH (DATE,GMT,SLAT,SLON,GLAT,GLON,AERR)
CALL GECS(GLAT,GLON,SLAT,SLON,ASAT,ZSAT)
CALL HOCS(DATE,GMT,GLAT,GLON,DEC,RA,ALT,AZI)
PRINT 200,DATE,TIME
200 FORMAT (//,5X,’DATE: ‘,A12,10X,’TIME:-‘,A10,'(GMT)’/)
PRINT 300
300 FORMAT (5X,’POSITION OF SUN :’/5X,17(‘=’)/)
PRINT 400,ECLON
400 FORMAT (5X,’ECLIPTIC LONGITUDE:’,F7.2)
PRINT 500,DEC,RA
500 FORMAT(5X,’DECLINATION:’F9.2/5X,’RIGHT ASCENTION:’,F9.2)
PRINT 600,ALT,AZI
600 FORMAT(5X,’ELEVATION:’,F9.2/5X,’AZIMUTH:’,F9.2/)
PRINT 700
700 FORMAT(5X,’POSITION OF THE MOON:’/5X,27(‘=’)/)
PRINT 800,SLAT,SLON
800 FORMAT(5X,’SUB-LUNAR POINT:(SLAT=’,F6.2,’,SLON=’,F7.2,’)’/)
PRINT 900, ASAT,ZSAT
900 FORMAT(5X,’ELEVATION:’,F9.2/5X,’AZIMUTH:’,F9.2/)
PRINT 905
905 FORMAT(5X,’POSITION OF THE SPECULAR POINT :’/5X,32(‘=’)/)
IF (AERR .EQ. .0D0) THEN
PRINT 910,GLAT,GLON
910 FORMAT (5X,’LATITUDE:’,F9.2/5X,’LONGITUDE:’,F7.2)
ELSE
PRINT 915
915 FORMAT(5X,’GLINT DOES NOT LIE WITH IN THE VIEW AREA’)
END IF
END

DOUBLE PRECISION FUNCTION HOURANG(DATE,GMT,GLON,RA)
IMPLICIT LOGICAL (A-Z)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CHARACTER*15 DATE
DOUBLE PRECISION GMT,GLON,RA,GLON1,HANG
DOUBLE PRECISION LST,GST
GLON1=GLON
IF (GLON1 .GT. .18d3) GLON1=GLON1-.36D3
lst=glon1/.15d2+gst(date,gmt)
IF (LST .GT. .24D2) LST=LST-.24D2
IF (LST .LT. .0D0) LST=LST+.24D2
HANG=(LST-RA/.15D2)
IF (HANG .LT. .0D0) HANG=HANG+.24D2
HOURANG=HANG
RETURN
END

SUBROUTINE NEWRAPH(DATE,GMT,SLAT,SLON,GLAT,GLON,AERR)
IMPLICIT logical (A-Z)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INTEGER KK
DOUBLE PRECISION HSAT,ERAD,PI,RAD,RADI,CDIFF,DEC,RA,SLAT,GLAT
DOUBLE PRECISION GMT,SANG,DIFF,SLON,GLON,GLAT1,AERR,GLON1,SLON1
DOUBLE PRECISION SLAT1,DEC1,DLAT,DLON,OLDGLAT1,jd
DOUBLE PRECISION OLDGLON1,GLON2,U1,U2,U3,U4,V1,V2,V3,HANG1,HANG
DOUBLE PRECISION F0,G0,DUDT1,DUDT2,DUDT3,DUDT4,DUDP1,DUDP2,DUDP3
DOUBLE PRECISION DUDP4,HOURANG,ATND,DVDT1,DVDT2,DVDT3,XX1,XX2
DOUBLE PRECISION DVDP1,DVDP2,DVDP3,T1,T2,T3,T4,TEMP1,TEMP2,TEMP3
DOUBLE PRECISION AK,AH,TEMP4,TEMP5,TEMP6,DFDT,DFDP,DGDT,DGDP
DOUBLE PRECISION ALATD,ALOND,ECLAT,ECLON
CHARACTER*15 DATE
PARAMETER (HSAT=.384d9,ERAD=.6378D7)
c PARAMETER (HSAT=.35786D8,ERAD=.6378D7)
PARAMETER (PI=.3141592654D1)
RAD=PI/.18D3
RADI=.18D3/PI
CDIFF=.1D1/.36D4*RAD
CALL ECCS(DATE,GMT,ECLAT,ECLON,jd)
CALL EQCS(DATE,GMT,ECLAT,ECLON,DEC,RA)
GLAT=(DEC+SLAT)/.2D1
SANG=(.12D2-GMT)*.15D2
IF (SANG.LT. .0D0) SANG=SANG+.36D3
DIFF=dabs(SANG-SLON)
IF (DIFF .GT. .18D3) THEN
DIFF=.36D3-DIFF
GLON=DMAX1(SANG,SLON)+DIFF/.2D1
ELSE
GLON=DMIN1(SANG,SLON)+DIFF/.2D1
END IF
IF (GLON .GT. .36D3) GLON=GLON-.36D3
IF (DIFF .GT. .12D3) THEN
AERR=-.1D1
RETURN
END IF
GLAT1=GLAT*RAD
GLON1=GLON*RAD
SLAT1=SLAT*RAD
SLON1=SLON*RAD
DEC1=DEC*RAD
KK=0
150 KK=KK+1
DLAT=GLAT1-SLAT1
DLON=GLON1-SLON1
OLDGLAT1=GLAT1
OLDGLON1=GLON1
GLON2=GLON1*RADI
HANG=HOURANG ( DATE , GMT , GLON2 , RA)
HANG1=HANG* .15D2*RAD
U1= ( HSAT+ERAD) *DCOS ( DLAT) *DCOS( DLON) -ERAD
U2=DSIN(GLAT1)*DSIN(DEC1)+DCOS(DEC1)*DCOS(GLAT1)*DCOS(HANG1)
U3=DSIN( DLON)
U4=-DCOS(DEC1)*DCOS(GLAT1)*DSIN(HANG1)
V1=HSAT**.2D1+.2D1*ERAD* (ERAD+HSAT)*(.1D1-DCOS(DLAT)*DCOS(DLON))
V2=DCOS ( DLON ) *DSIN (DLAT)
V3=DSIN(DEC1)-DSIN(GLAT1)*U2
F0=DASIN(U1/SQRT(V1) )-DASIN(U2)
G0=ATND(U3 ,V2 ) -ATND(U4 ,V3)
DUDT1=- ( ERAD+HSAT) *DSIN ( DLAT ) *DCOS( DLON)
DUDT2=DCOS(GLAT1)*DSIN(DEC1)-DSIN(GLAT1)*DCOS(DEC1)*DCOS(HANG1)
DUDT3= .0D0
DUDT4=DCOS ( DEC1 ) *DSIN(GLAT1 ) *DSIN(HANG1)
DUDP1=- ( ERAD+HSAT) *DCOS ( DLAT ) *DSIN ( DLON)
DUDP2= .0D0
DUDP3=DCOS( DLON)
DUDP4= .0d0
DVDT1= .2D1*ERAD* (HSAT+ERAD) *DSIN ( DLAT) *DCOS (DLON)
DVDT2=DCOS( DLAT) *DCOS ( DLON)
XX1=-.2D1*DSIN(GLAT1 )*DCOS(GLAT1)*DSIN(DEC1 )
XX2=DCOS(GLAT1) *DSIn(GLAT1)*DCOS(DEC1)*DCOS(HANG1)
DVDT3=XX1+XX2
DVDP1=.2D1*ERAD* (HSAT+ERAD) *DCOS( DLAT) *DSIN (DLON )
DVDP2=-DSIN( DLAT) *DSIN ( DLON)
DVDP3=.0D0
T1=.1D1/(SQRT(.1D1-(U1/SQRT(V1))**.2D1))
T2=.1D1/(SQRT(.1D1-U2**.2D1))
T3=.1D1/(.1D1+(U3/V2)**.2D1)
T4=.1D1/(.1D1+(U4/V3)**.2D1)
TEMP1=(SQRT(V1)*DUDT1-U1*DVDT1/(SQRT(V1)*.2D1))/V1
TEMP2=(SQRT(V1)*DUDP1-U1*DVDP1/(SQRT(V1)*.2D1))/V1
TEMP3=(V2*DUDT3-U3*DVDT2)/V2**.2D1
TEMP4=(V2*DUDP3-U3*DVDP2)/V2**.2D1
TEMP5=(V3*DUDT4-U4*DVDT3)/V3**.2D1
TEMP6=(V3*DUDP4-U4*DVDP3)/V3**.2D1
DFDT=T1*TEMP1-T2*DUDT2
DFDP=T1*TEMP2-T2*DUDP2
DGDT=T3*TEMP3+T4*TEMP5*DVDT3
DGDP=T3*TEMP4+T4*TEMP6*DVDP3
AK=(F0*DGDT-G0*DFDT)/(DGDP*DFDT-DFDP*DGDT)
AH=-(AK*DFDP+F0)/DFDT
GLAT1=GLAT1+AH
GLON1=GLON1+AK
IF (GLON1.LT. .0D0) GLON1=GLON1+.2D1*PI
IF (GLON1.GT. .2D1*PI) GLON1=GLON1-.2D1*PI
IF (GLAT1.LT. -PI/.2D1) GLAT1=-PI-GLAT1
IF (GLAT1.GT. PI/.2D1) GLAT1=PI-GLAT1
ALATD=dabs(GLAT1-OLDGLAT1)
ALOND=dabs(GLON1-OLDGLON1)
IF (ALOND .GT.PI) ALOND=.2D1*PI-ALOND
IF (KK .GT.150) THEN
AERR=-.2D0
RETURN
END IF
IF (ALOND .GT.CDIFF .AND.ALATD .GT.CDIFF) GOTO 150
GLAT=GLAT1*RADI
GLON=GLON1*RADI
AERR=.0D0
RETURN
END

SUBROUTINE ECCS (DATE , GMT , ECLAT , ECLON, JD)
IMPLICIT LOGICAL (A-Z)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION GMT , ECLAT , ECLON, JD,ELE,ELP,ECC,PI
DOUBLE PRECISION EC,JDEP,RAD,AJDATE,DAYS,AN,M,AM
CHARACTER*15 DATE
PARAMETER (ELE=.279403303D3 ,ELP=.282768422D3)
PARAMETER (ECC=.167131D-1 ,PI=.3141592654D1)
RAD=PI/.18D3
JDEP=.24478915D7
JD=AJDATE(DATE,GMT)
DAYS=JD-JDEP
AN=(.36D3/.365242191D3)*DAYS
M=DINT(AN/.36D3)
AN=AN-(M*.36D3)
AM=AN+ELE-ELP
IF (AM .LT. .0D0) AM=AM+.36D3
EC=( .36D3 /PI ) *ECC*DSIN(Am*RAD)
ECLON=AN+EC+ELE
IF (ECLON .GT. .36d3) ECLON=ECLON-.36D3
IF (ECLON .LT. .0D0) ECLON=ECLON+.36D3
ECLAT=.0D0
RETURN
END

SUBROUTINE EQCS ( DATE, GMT , ECLAT , ECLON , DEC,RA)
IMPLICIT LOGICAL (A-Z)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CHARACTER*15 DATE
DOUBLE PRECISION JD,JDEP,PI,GMT , ECLAT , ECLON , DEC,RA
DOUBLE PRECISION RAD,RADI,T,EP,ECLON1,ECLAT1,TEMP1,TEMP2
DOUBLE PRECISION AJDATE,DEP,ANR,DNR,ATND
PARAMETER (PI=.3141592654D1)
RAD=PI/.18D3
RADI=.18D3/PI
JDEP=.2451545D7
JD=AJDATE(DATE,GMT)
T=(JD-JDEP)/.36525D5
DEP=.46815D2*T+.6D-3*T**.2D1-.181D-2*T**.3D1
EP=.23439292D2-DEP/.36D4
ECLON1=ECLON*RAD
ECLAT1=ECLAT*RAD
EP=EP*RAD
TEMP1=DSIN(ECLAT1)*DCOS(EP)
TEMP2=DCOS(ECLAT1)*DSIN(EP)*DSIN(ECLON1)
DEC=DASIN(TEMP1+TEMP2)*RADI
ANR=DSIN(ECLON1)*DCOS(EP)-DTAN(ECLAT1)*DSIN(EP)
DNR=DCOS(ECLON1)
RA=ATND(ANR,DNR)*RADI
RETURN
END

SUBROUTINE HOCS(DATE,GMT,GLAT,GLON,DEC,RA,ALT,AZI)
IMPLICIT LOGICAL (A-Z)
c IMPLICIT DOUBLE PRECISION (A-Z)
DOUBLE PRECISION GMT,GLAT,GLON,DEC,RA,ALT,AZI,U1,U2,V2
DOUBLE PRECISION PI,RAD,RADI,HANG,HANG1,GLAT1,DEC1
DOUBLE PRECISION HOURANG,ATND
PARAMETER (PI=.3141592654D1)
CHARACTER*15 DATE
RAD=PI/.18D3
RADI=.18D3/PI
HANG=HOURANG(DATE,GMT,GLON,RA)
HANG1=HANG*.15D2*RAD
DEC1=DEC*RAD
GLAT1=GLAT*RAD
U1=DSIN(DEC1)*DSIN(GLAT1)+DCOS(DEC1)*DCOS(GLAT1)*DCOS(HANG1)
ALT=DASIN(U1)*RADI
U2=-DCOS(DEC1)*DSIN(HANG1)*DCOS(GLAT1)
V2=DSIN(DEC1)-DSIN(GLAT1)*DSIN(ALT*RAD)
AZI=ATND(U2,V2)*RADI
RETURN
END

DOUBLE PRECISION FUNCTION GST(DATE,GMT)
IMPLICIT LOGICAL (A-Z)
c IMPLICIT DOUBLE PRECISION (A-Z)
CHARACTER*15 DATE
DOUBLE PRECISION JD,JDEP,T,T0,TT,GST1,GMT,AJDATE
JDEP=.2451545D7
JD=AJDATE(DATE,.0D0)
T=(JD-JDEP)/.36525D5
T0=.6697374558D1+.2400051336D4*T+.25862D-4*T**.2D1
TT=DINT(abs(T0)/.24D2)+.1D1
T0=T0+TT*.24D2
GST1=GMT*.1002737909D1+T0
TT=DINT(GST1/.24D2)
GST1=GST1-TT*.24D2
GST=GST1
RETURN
END

DOUBLE PRECISION FUNCTION ATND(Y,X)
IMPLICIT LOGICAL (A-Z)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION x,y,pi,dir
PARAMETER (PI=.3141592654D1)
IF (X .EQ. .0D0 .OR.Y .EQ. .0D0) THEN
IF (X .EQ. .0D0 .AND.Y .EQ. .0D0) DIR=.0D0
IF (X .EQ. .0D0 .AND.Y .GT. .0D0) DIR=PI/.2d1
IF (X .EQ. .0D0 .AND.Y .LT. .0D0) DIR=.15D1*PI
IF (Y .EQ. .0D0 .AND.X .GT. .0D0) DIR=.0D0
IF (Y .EQ. .0D0 .AND.X .LT. .0D0) DIR=PI
ELSE
DIR=ATAN(Y/X)
IF (Y .GT. .0D0 .AND.X .LT. .0D0 ) DIR=DIR+PI
IF (Y .LT. .0D0 .AND.X .LT. .0D0 ) DIR=DIR-PI
IF (DIR .LT. .0D0) DIR=DIR+.2D1*PI
END IF
ATND=DIR
RETURN
END

SUBROUTINE GECS(GLAT,GLON,SLAT,SLON,ASAT,ZSAT)
IMPLICIT LOGICAL (A-Z)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION HSAT,ERAD,PI,GLAT,GLON,SLAT,SLON,ASAT,ZSAT
DOUBLE PRECISION RAD,RADI,GLAT1,GLON1,SLAT1,SLON1,DLAT,DLON
DOUBLE PRECISION U1,V1,U2,V2,ATND
PARAMETER (HSAT=.384d9,ERAD=.6378D7)
c PARAMETER (HSAT=.35786D8,ERAD=.6378D7)
PARAMETER (PI=.3141592654D1)
RAD=PI/.18D3
RADI=.18D3/PI
GLAT1=GLAT*RAD
GLON1=GLON*RAD
SLAT1=SLAT*RAD
SLON1=SLON*RAD
DLAT=(GLAT1-SLAT1)
DLON=(GLON1-SLON1)
U1=(HSAT+ERAD)*DCOS(DLAT)*DCOS(DLON)-ERAD
V1=HSAT**.2D1+.2D1*ERAD*(ERAD+HSAT)*(.1D1-DCOS(DLAT)*DCOS(DLON))
U2=DSIN(DLON)
V2=DCOS(DLON)*DSIN(DLAT)
ASAT=DASIN(U1/SQRT(V1))*RADI
ZSAT=(ATND(U2,V2)+PI)*RADI
IF (ZSAT .GT. .36D3) ZSAT=ZSAT-.36D3
RETURN
END

DOUBLE PRECISION FUNCTION AJDATE(DATE,GMT)
IMPLICIT LOGICAL (A-Z)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION A,B,C,D,MONTH,JDBC,GMT,DAY,YEAR
PARAMETER (JDBC=.17209945D7)
CHARACTER*15 DATE,INSTR
WRITE(INSTR,100) DATE
100 FORMAT (A15)
READ (INSTR,200) DAY,MONTH,YEAR
200 FORMAT (2(F2.0,1X),F4.0)
DAY=DAY+GMT/.24D2
IF (MONTH .LT. .3D1) THEN
YEAR=YEAR-.1D1
MONTH=MONTH+.12D2
END IF
A=DINT(YEAR/.1D3)
B=.2D1-A+DINT(A/.4D1)
C=DINT(.36525D3*YEAR)
D=DINT(.306001D2*(MONTH+.1D1))
ajdate=b+c+d+day+jdbc
return
end

I have made an attempt to figure out what the noise levels are at various fluxes in individual frames.

28 images of the nearly-full moon were used. They were very carefully aligned (using subpixel rebinning and the iraf task imshift). The shifts were determined with a fortran program which centered the moon using center of light — this was written because no iraf task I could find could easily align such images without reference stars.

Below I show in green circles the regions that were used to measure the standard deviation of pixels in a 10×10 box on the lunar face. These were chosen because they are relatively flat and have a range of luminosities. All the pixels in the scattered light off the lunar edge were also used, in a sector starting from the center between the position angles 40 and 50 degrees (i.e. off to the upper right at 45 degrees).

The standard deviation of the selected pixel positions over all the 28 exposures was computed, and a plot made of the sqrt of the mean counts in the pixel versus the standard deviation of the counts.

If the noise is Poissonian, then these are related 1:1. The plot below shows that this is only approximately the case (solid line). A closer model is that the noise is about half of what is expected (dotted line). Why the noise is sub-Poissonian is open to speculation, but must be due to some slight smoothing we don’t account for (possibly iraf introduces some smoothing when rebinning the data for the subpixel shifting).

The points at the lower left are those from the scattered light off the lunar edge. The other points are the pixels in the 5 small patches indicated by the green circles in the previous figure.

It’s difficult to fill in more points because it’s hard to find flat regions of the appropriate luminosity. Flat regions are needed because the technique is very sensitive to slight misalignments in the images (if one examines regions where the luminosity is changing rapidly — such as near the lunar edge, craters, strirations etc on the lunar face — the standard deviation goes haywire).

Assuming Poissoninan noise in test data generated from synthetic moons appears to be a conservative option, as it is pretty close to the correct value and possibly an overestimate, so currently I’d recommend that.

Update and error correction: The ADU conversion for this chip is approx. 3.8. The statistics above are were done in ADU, not photons. This of course accounts for the factor of two (=sqrt(3.8)) problem above. Thanks to Peter for spotting this error.

Using photons instead of electrons gives a pretty close fit t the 1:1 relationship. We could try this for a range of lunar exposure times, using the same flat regions on the moon and the halo off the edge, to fill in the gaps for a range of luminosities.

Getting the setting Moon in clear skies.

Plenty of ‘filter stick’ and it comes in sets – when we change filters the shutter ‘unsticks’ and we are good for one set of images, then we change filter and there is the chance that we get good data or bad data.

We also see a periodic shifting of the Moon inside the frame. Either the JPL table of Moon positions is slightly bad, or the interpolation is done.

Question: Is the table of Moon positions interpolated or do we go for ‘nearest neighbour’?

See the telescope and the Moon!