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$ AbstractCC This subroutine is called by LITHOS and LITHOSP if ALPAD=.TRUE. inC 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.CC The input weighting factor SK = SCALE**2/(SCALE**2+PICRES**2) approachesC 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).CC 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.CC$ DisclaimerCCC$ Required_ReadingCC R.W. Gaskell, et.al, "Characterizing and navigating small bodiesC with imaging data", Meteoritics & Planetary Science 43,C Nr 6, 1049-1061 (2008)CC$ 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/OCC Variable I/O DescriptionC -------- --- --------------------------------------------------C PMX I Maximum number of images contributing to a mapletC NTMP I Maximum maplet pixel rangeC 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 gridC DN1 I Data grid with extracted and predicted maplet brightnessC SP I Sun vector for each image in maplet frameC CP I Camera vector for each image in maplet frameC PUSE I .TRUE. if image used in maplet constructionC LAMBDA I Current brightness coefficient for each imageC PHI I Current background brightness for each picture C SK I A weighting factor emphasizing higher res imagesC TUSE I .TRUE. means maplet data exists at that pixelC TMPL I/O Slope and albedo deviation at each maplet pixelCC$ File_I/OCC Filename I/O DescriptionC ---------------------------- --- -------------------------------C NoneCC$ RestrictionsCCC$ Software_DocumentationCC OSIRIS-REx Stereophotoclinometry Software Design DocumentC OSIRIS-REx Stereophotoclinometry Software User's GuideCC$ Author_and_InstitutionCC R.W. Gaskell (PSI)CC$ VersionCCCC$ SPC_functions_calledC ILLUMCC$ SPC_subroutines_calledC NoneCC$ SPICELIB_functions_calledC RPD % Conversion for radians/degreesC VDOTCC$ SPICELIB_subroutines_calledC NoneC$ Called_by_SPC_ProgramsC LITHOSC 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