== 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

{{attachment:coverage_p.jpg}}



{{{
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



}}}