v3.0.3 - 11 Mar 2016

C$ Procedure SUBROUTINE FIND_LAMBDA_PIC(NTMP,QSZ,DNX, SP,CP,TUSE,TMPL, . LAMBDA,PHI,CHI,CHI0) C$ Abstract C C The predicted brightness for a pixel in a maplet is given by C C B = LAMBDA*ETA*ILLUM(COS(I),COS(E),PHASE)+PHI C C where I and E are incidence and emission angles, ETA is relative C albedo, normalized to 1 across the maplet. Since no attempt is made C to account for absolute radiometry, there is a constant lambda that C serves to relate the topography manifested by the slopes leading to C variations in I and E to the brightness measured by the imaging data. C Clearly, a smaller value of lambda requires a more robust topography C than a larger one, so we must have some initial idea of the topography C to determine lambda. This is accomplished by making large, low C resolution maplets to start, that owe most of their relief to the C curvature of the body from ordinary stereography. Once lambda is C determined, it sets the scale for the higher frequency topography in C the maplet. This topography, in turn, sets the scale for the higher C resolution maplets that overlap the original maplet and their even C higher frequency topography. The topography scale of any maplet is C constantly refined by contributions such as differential stereo, limb C contributions and randomly selected heights from overlapping maplets. C C This subroutine estimates the value oflambda in the brightness function C from the observed brightness distribution and the nominal topography. It C can also solve for a constant background in the brightness such as was C observed due to atmospheric scattering in Viking Orbiter images of Mars. C This latter PHI term is usually ignored by setting the parameter CHI0 = 1 C in INIT_LITHOS.TXT. 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 NTMP DOUBLE PRECISION SP(3) DOUBLE PRECISION CP(3) DOUBLE PRECISION ALPHA DOUBLE PRECISION RPD DOUBLE PRECISION VDOT DOUBLE PRECISION ILLUM DOUBLE PRECISION GAMMA DOUBLE PRECISION ETA DOUBLE PRECISION LAMBDA DOUBLE PRECISION PHI DOUBLE PRECISION CHI DOUBLE PRECISION CHI0 DOUBLE PRECISION X1 DOUBLE PRECISION X2 DOUBLE PRECISION Z0 DOUBLE PRECISION Z1 DOUBLE PRECISION Z2 DOUBLE PRECISION ZI DOUBLE PRECISION ZB DOUBLE PRECISION ZBI DOUBLE PRECISION ZBB DOUBLE PRECISION ZII REAL*4 DNX(-NTMP:NTMP,-NTMP:NTMP) REAL*4 TMPL(-NTMP:NTMP,-NTMP:NTMP,3) INTEGER QSZ INTEGER I INTEGER J INTEGER N1 LOGICAL TUSE(-NTMP:NTMP,-NTMP:NTMP) C$ Variable_I/O C C Variable I/O Description C -------- --- -------------------------------------------------- C NTMP C QSZ C DNX C SP C CP C TUSE C TMPL C LAMBDA C PHI C CHI C CHI0 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 C VDOT C C$ SPICELIB_subroutines_called C None C C$ Called_by_SPC_Programs C AUTOREGISTER C SUBROUTINE FIND_LAMBDA LAMBDA=1 PHI=0 Z0=1 ZI=0 ZB=0 ZBI=0 ZBB=0 ZII=0 N1=0 % Do per pixel DO I=-QSZ,QSZ DO J=-QSZ,QSZ IF((DNX(I,J).GT.0).AND.TUSE(I,J)) THEN X1=DNX(I,J) Z0=MIN(X1,Z0) GAMMA=SQRT(1+TMPL(I,J,1)**2+TMPL(I,J,2)**2) Z1=(SP(3)+TMPL(I,J,1)*SP(1)+TMPL(I,J,2)*SP(2))/GAMMA Z2=(CP(3)+TMPL(I,J,1)*CP(1)+TMPL(I,J,2)*CP(2))/GAMMA ETA=1+TMPL(I,J,3) ALPHA=ACOS(VDOT(CP,SP))/RPD() % phase X2=ETA*ILLUM(Z1,Z2,ALPHA) ZI=ZI+X1 ZB=ZB+X2 ZBI=ZBI+X1*X2 ZBB=ZBB+X2**2 ZII=ZII+X1**2 N1=N1+1 ENDIF ENDDO ENDDO IF((N1.GT.0).AND.(ZB.GT.0)) THEN ZI=ZI/N1 ZB=ZB/N1 ZBI=ZBI/N1 ZBB=ZBB/N1 ZII=ZII/N1 LAMBDA=ZI/ZB PHI=0 CHI=0 IF((ZBB-ZB**2.GT.0.001).AND.(ZII-ZI**2.GT.0.001)) THEN CHI=(ZBI-ZB*ZI)/SQRT((ZBB-ZB**2)*(ZII-ZI**2)) IF(CHI.GT.CHI0) THEN Z1=(ZI*ZBB-ZB*ZBI)/(ZBB-ZB**2) IF(Z1.GT.0) THEN PHI=MIN(Z0,Z1) LAMBDA=(ZI-PHI)/ZB ENDIF ENDIF ENDIF ENDIF RETURN END