find_albedo here
From v3.0.3 - May 11 2016
WGT=SK(K)*(1-MU*SP(3,K)**2)/LAMBDA(K)**2 K is the image Mu = 0.75 SP is sun position Lambda(K) Current brightness coefficient for each image
The reduction base on incidence angle
C$ Procedure
SUBROUTINE FIND_ALBEDO(PMX,NTMP,NPIX,QSZ,NSZ,DN1,SP,CP,
. PUSE,LAMBDA,PHI,SK,TUSE,TMPL)
C$ Abstract
C
C This subroutine is called by LITHOS and LITHOSP if ALPAD=.TRUE. in
C INIT_LITHOS.TXT. It is intended to give a better estimate for the
C albedo at each pixel of the maplet by more appropriate weighting of
C the data.
C
C The input weighting factor SK = SCALE**2/(SCALE**2+PICRES**2) approaches
C 1 for very high resolution images and 0 for very low resolution images.
C It is then multiplied by a factor [1-0.75cos(i)] that de-emphasizes high
C local sun angles and by B^4/(B^4+nu^4) that de-emphasizes brightness less
C than nu=.05 (5% of full scale).
C
C The sum of the weighted squared residuals between observed and predicted
C brightness for each image contributing to each pixel of the maplet are
C minimized to obtain a correction to the albedo.
C
C$ Disclaimer
C
C
C$ Required_Reading
C
C R.W. Gaskell, et.al, "Characterizing and navigating small bodies
C with imaging data", Meteoritics & Planetary Science 43,
C Nr 6, 1049-1061 (2008)
C
C$ Declarations
IMPLICIT NONE
INTEGER PMX % Max number of images that SPC can hold
INTEGER NTMP % Variable for the size of maplet
DOUBLE PRECISION SP(3,PMX)
DOUBLE PRECISION CP(3,PMX)
DOUBLE PRECISION ALPHA
DOUBLE PRECISION RPD
DOUBLE PRECISION VDOT
DOUBLE PRECISION ILLUM
DOUBLE PRECISION GAMMA
DOUBLE PRECISION ETA
DOUBLE PRECISION LAMBDA(PMX)
DOUBLE PRECISION PHI(PMX)
DOUBLE PRECISION SK(PMX)
DOUBLE PRECISION RSD(PMX)
DOUBLE PRECISION Z1, Z2, Z3, Z5, Z6
DOUBLE PRECISION W, M, P, WGT, MU, NU
REAL*4 DN1(5000,5000)
REAL*4 TMPL(-NTMP:NTMP,-NTMP:NTMP,3)
INTEGER NPIX % Number of images contributing to maplet
INTEGER QSZ % Maplet pixel range from -qsz to qsz
INTEGER NSZ % Number of maplets in one row of data grid
INTEGER I
INTEGER J
INTEGER K
INTEGER N
INTEGER K1
INTEGER N1
INTEGER N2
INTEGER I0
INTEGER J0
LOGICAL TUSE(-NTMP:NTMP,-NTMP:NTMP)
LOGICAL PUSE(PMX)
C$ Variable_I/O
C
C Variable I/O Description
C -------- --- --------------------------------------------------
C PMX I Maximum number of images contributing to a maplet
C NTMP I Maximum maplet pixel range
C NPIX I Number of images contributing to maplet.
C QSZ I Maplet pixel range from -qsz to qsz.
C NSZ I Number of maplets in one row of data grid
C DN1 I Data grid with extracted and predicted maplet brightness
C SP I Sun vector for each image in maplet frame
C CP I Camera vector for each image in maplet frame
C PUSE I .TRUE. if image used in maplet construction
C LAMBDA I Current brightness coefficient for each image
C PHI I Current background brightness for each picture
C SK I A weighting factor emphasizing higher res images
C TUSE I .TRUE. means maplet data exists at that pixel
C TMPL I/O Slope and albedo deviation at each maplet pixel
C
C$ File_I/O
C
C Filename I/O Description
C ---------------------------- --- -------------------------------
C None
C
C$ Restrictions
C
C
C$ Software_Documentation
C
C OSIRIS-REx Stereophotoclinometry Software Design Document
C OSIRIS-REx Stereophotoclinometry Software User's Guide
C
C$ Author_and_Institution
C
C R.W. Gaskell (PSI)
C
C$ Version
C
C
C
C$ SPC_functions_called
C ILLUM
C
C$ SPC_subroutines_called
C None
C
C$ SPICELIB_functions_called
C RPD % Conversion for radians/degrees
C VDOT
C
C$ SPICELIB_subroutines_called
C None
C$ Called_by_SPC_Programs
C LITHOS
C
MU=0.75
NU=0.05
DO I=-QSZ,QSZ
DO J=-QSZ,QSZ
% Assuming TMPL (template) with 1 and 2 being slope x/y and 3 being albedo
% Assuming gamma is slope
% Assuming eta is albedo (why do we add 1?)
IF(TUSE(I,J)) THEN
GAMMA=SQRT(1+TMPL(I,J,1)**2+TMPL(I,J,2)**2)
ETA=1+TMPL(I,J,3)
M=0
W=0
N=0
N1=0
% For each image
DO K=1,NPIX
IF(PUSE(K)) THEN
N=N+1 % counts the number of images used
K1=(N-1)/NSZ % scaled for the number of maplets
N2=N-1-NSZ*K1
I0=(QSZ+1)*(1+2*N2)
J0=(QSZ+1)*(1+4*K1)
Z5=(SP(3,K)+TMPL(I,J,1)*SP(1,K)+TMPL(I,J,2)*SP(2,K))/GAMMA
Z6=(CP(3,K)+TMPL(I,J,1)*CP(1,K)+TMPL(I,J,2)*CP(2,K))/GAMMA
% If valid
IF((Z5.GT.0).AND.(Z6.GT.0)) THEN
IF(DN1(I+I0,J+J0).GT.0) THEN
% Alpha is phase?
ALPHA=ACOS(VDOT(CP(1,K),SP(1,K)))/RPD()
% RSD is residual? For the whole image? It gets overwritten
RSD(K)=DN1(I+I0,J+J0)
. -LAMBDA(K)*ETA*ILLUM(Z5,Z6,ALPHA)-PHI(K)
P=LAMBDA(K)*ILLUM(Z5,Z6,ALPHA)
WGT=SK(K)*(1-MU*SP(3,K)**2)/LAMBDA(K)**2
Z2=DN1(I+I0,J+J0)
WGT=WGT*Z2**4/(Z2**4+NU**4)
% What is M and W?
M=M+P*P*WGT
W=W+RSD(K)*P*WGT
N1=N1+1
ENDIF
ENDIF
ENDIF
ENDDO
IF(N1.GE.1) TMPL(I,J,3)=TMPL(I,J,3)+SNGL(W/M)
ENDIF
ENDDO
ENDDO
Z1= 1.D10
Z2=-1.D10
Z3=0
N1=0
DO I=-QSZ,QSZ
DO J=-QSZ,QSZ
IF(TUSE(I,J)) THEN
Z1=MIN(Z1,TMPL(I,J,3))
Z2=MAX(Z2,TMPL(I,J,3))
Z3=Z3+TMPL(I,J,3)
N1=N1+1
ENDIF
ENDDO
ENDDO
Z3=Z3/N1
Z1=1.01*MAX(1.D0,(Z2-Z3),(Z3-Z1))
DO I=-QSZ,QSZ
DO J=-QSZ,QSZ
IF(TUSE(I,J)) THEN
TMPL(I,J,3)=(TMPL(I,J,3)-SNGL(Z3))/SNGL(Z1)
ENDIF
ENDDO
ENDDO
RETURN
END