$ cat earthshine.m 
function alb = earthshine(refl, degrees);
sunWidth = pi/4
sunAngle = -pi/4
satShift = pi - sunWidth
f = degrees * pi/180 - pi / 2 + sunWidth / 2
earthR = 6371 * 1000
satH = 400*1000*1000 + earthR
sunH = 150*1000*1000*1000
sat = [cos(f+satShift)*satH,sin(f+satShift)*satH,0]
sunX = sunH
sunZ = 0
sunX = sunX * cos(sunAngle) - sunZ * sin(sunAngle)
sunZ = sunX * sin(sunAngle) + sunZ * cos(sunAngle)
sun = [cos(f)*sunX,sin(f)*sunX,sunZ]
a = albedo(sat, sun, refl) % 's' 'p') 
alb = sum(sum(a))
return

filesJune = [
, "ga050501.mat"
, "ga050502.mat"
, "ga050503.mat"
, "ga050504.mat"
, "ga050505.mat"
, "ga050506.mat"
, "ga050507.mat"
, "ga050508.mat"
, "ga050509.mat"
, "ga050510.mat"
, "ga050511.mat"
, "ga050512.mat"
, "ga050513.mat"
, "ga050514.mat"
, "ga050515.mat"
, "ga050516.mat"
, "ga050517.mat"
, "ga050518.mat"
, "ga050519.mat"
, "ga050520.mat"
, "ga050521.mat"
, "ga050522.mat"
, "ga050523.mat"
, "ga050524.mat"
, "ga050525.mat"
, "ga050526.mat"
, "ga050527.mat"
, "ga050528.mat"
, "ga050529.mat"
, "ga050530.mat"
, "ga050531.mat"
]

filesDec = [
, "ga051201.mat"
, "ga051202.mat"
, "ga051203.mat"
, "ga051204.mat"
, "ga051205.mat"
, "ga051206.mat"
, "ga051207.mat"
, "ga051208.mat"
, "ga051209.mat"
, "ga051210.mat"
, "ga051211.mat"
, "ga051212.mat"
, "ga051213.mat"
, "ga051214.mat"
, "ga051215.mat"
, "ga051216.mat"
, "ga051217.mat"
, "ga051218.mat"
, "ga051219.mat"
, "ga051220.mat"
, "ga051221.mat"
, "ga051222.mat"
, "ga051223.mat"
, "ga051224.mat"
, "ga051225.mat"
, "ga051226.mat"
, "ga051227.mat"
, "ga051228.mat"
, "ga051229.mat"
, "ga051230.mat"
, "ga051231.mat"
]

reflJune = avg_data(filesJune)
reflDec  = avg_data(filesDec)
june = avg_alb(reflJune) % set earthshine's sunAngle = +pi/4
dec = avg_alb(reflDec)   % set earthshine's sunAngle = -pi/4

plot(circshift(june,180))
plot(circshift(dec,180))

degrees = 45
sunWidth = pi/4
sunAngle = -pi/4
satShift = pi - sunWidth
f = degrees * pi/180 - pi / 2 + sunWidth / 2
earthR = 6371 * 1000
satH = 400*1000*1000 + earthR
sunH = 150*1000*1000*1000
sat = [cos(f+satShift)*satH,sin(f+satShift)*satH,0]
sunX = sunH
sunZ = 0
sunX = sunX * cos(sunAngle) - sunZ * sin(sunAngle)
sunZ = sunX * sin(sunAngle) + sunZ * cos(sunAngle)
sun = [cos(f)*sunX,sin(f)*sunX,sunZ]
a = albedo(sat, sun, reflDec, 'p')