v 3.0.3 - May 11 2016

C$ Procedure SUBROUTINE FIND_TMPL(PMX,NPIX,RSD,SP,CP,APR,LAMBDA,WGT, . T1,T2,T3,USE,PUSE,TFND,RSDSQ,WSM,S2S)C$ AbstractCC This subroutine performs the estimation of slope and relative albedo C for each pont in the maplet, using residuals and weights from the C FIND_SLOPES subroutine. With T1 and T2 the slopes ant 1+T3 the relativeC albedo, we are minimizing:CC ∑W(k)(O - P - dT1∂P/∂T1 - dT2∂P/∂T2 - dT3∂P/∂T3)^2C kCC where the sum is over images, O is the observed brightness at a pixel andC P is the predicted brightness. The weights W(k) and residuals O-P areC carried in from the FIND_SLOPES subroutine. The partials with respect to C T1 and T2 are determined numerically.CC There is no reason why these two subroutines could not be combined, butC for the moment we are taking an "If it ain't broke ..." attitude. 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 DOUBLE PRECISION MAT(6,6) DOUBLE PRECISION IMAT(6,6) DOUBLE PRECISION SMAT(6,6) DOUBLE PRECISION W(3) DOUBLE PRECISION PAR(3) DOUBLE PRECISION SP(3,PMX) DOUBLE PRECISION CP(3,PMX) DOUBLE PRECISION APR(3) DOUBLE PRECISION ALPHA % phase angle DOUBLE PRECISION RPD DOUBLE PRECISION VDOT DOUBLE PRECISION RSD(PMX) DOUBLE PRECISION LAMBDA(PMX) DOUBLE PRECISION WGT(PMX) DOUBLE PRECISION ILLUM DOUBLE PRECISION GAMMA DOUBLE PRECISION ETA DOUBLE PRECISION WSM DOUBLE PRECISION Z2, Z3, Z4, Z5, Z6 DOUBLE PRECISION S2S DOUBLE PRECISION MXSLP DOUBLE PRECISION RSDSQ REAL*4 T1 REAL*4 T2 REAL*4 T3 INTEGER I1 INTEGER J1 INTEGER K INTEGER N1 INTEGER NPIX CHARACTER*80 LINE LOGICAL USE(PMX) LOGICAL PUSE(PMX) LOGICAL TFND LOGICAL EX DOUBLE PRECISION DELTA PARAMETER (DELTA=1.D-6) INTEGER*4 IFF SAVE IFF, MXSLP DATA IFF/0/ C$ Variable_I/OCC Variable I/O DescriptionC -------- --- --------------------------------------------------C PMXC NPIXC RSDC SPC CPC APRC LAMBDAC WGTC T1C T2C T3C USEC PUSEC TFNDC RSDSQC WSMC S2SCC$ File_I/OCC Filename I/O DescriptionC ---------------------------- --- -------------------------------C INIT_LITHOS.TXT I Setup file for the SPC toolkit CC$ 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 NoneCC$ SPC_subroutines_calledC INVERTNCC$ SPICELIB_functions_calledC NoneCC$ SPICELIB_subroutines_calledC NoneCC$ Called_by_SPC_ProgramsC LITHOSC SUBROUTINE FIND_SLOPESC % This is only done once IF(IFF.EQ.0) THEN IFF=1 MXSLP=0 INQUIRE(FILE='INIT_LITHOS.TXT',EXIST=EX) IF(EX) THEN OPEN(UNIT=25,FILE='INIT_LITHOS.TXT',STATUS='OLD')13 CONTINUE READ(25,FMT='(A80)') LINE IF(LINE(1:3).NE.'END') THEN IF(LINE(1:6).EQ.'MXSLP=') READ(LINE(7:80),*) MXSLP GO TO 13 ENDIF CLOSE(UNIT=25) ENDIF ENDIF DO I1=1,3 DO J1=1,3 MAT(I1,J1)=0 ENDDO W(I1)=0 ENDDO RSDSQ=0 WSM=0 N1=0 GAMMA=SQRT(1+T1**2+T2**2) ETA=1+T3 DO K=1,NPIX IF(USE(K).AND.PUSE(K)) THEN Z5=(SP(3,K)+T1*SP(1,K)+T2*SP(2,K))/GAMMA Z6=(CP(3,K)+T1*CP(1,K)+T2*CP(2,K))/GAMMA IF((Z5.GT.0).AND.(Z6.GT.0)) THEN ALPHA=ACOS(VDOT(CP(1,K),SP(1,K)))/RPD() % phase Z3=(ILLUM(Z5+DELTA,Z6,ALPHA)-ILLUM(Z5-DELTA,Z6,ALPHA)) . /(2*DELTA) Z4=(ILLUM(Z5,Z6+DELTA,ALPHA)-ILLUM(Z5,Z6-DELTA,ALPHA)) . /(2*DELTA) PAR(1)=(SP(1,K)/GAMMA-Z5*T1/GAMMA**2)*Z3 . +(CP(1,K)/GAMMA-Z6*T1/GAMMA**2)*Z4 PAR(2)=(SP(2,K)/GAMMA-Z5*T2/GAMMA**2)*Z3 . +(CP(2,K)/GAMMA-Z6*T2/GAMMA**2)*Z4 PAR(1)=LAMBDA(K)*ETA*PAR(1) PAR(2)=LAMBDA(K)*ETA*PAR(2) PAR(3)=LAMBDA(K)*ILLUM(Z5,Z6,ALPHA) DO I1=1,3 DO J1=1,3 MAT(I1,J1)=MAT(I1,J1)+PAR(I1)*PAR(J1)*WGT(K) ENDDO ENDDO DO I1=1,3 W(I1)=W(I1)+RSD(K)*PAR(I1)*WGT(K) ENDDO RSDSQ=RSDSQ+RSD(K)**2*WGT(K) WSM=WSM+WGT(K) N1=N1+1 ENDIF ENDIF ENDDO S2S=0 IF(N1.GT.2) THEN SMAT(1,1)=MAT(1,1)+0.01 SMAT(1,2)=MAT(1,2) SMAT(2,1)=MAT(2,1) SMAT(2,2)=MAT(2,2)+0.01 CALL INVERTN(2,SMAT,IMAT) S2S=(IMAT(1,1)+IMAT(2,2))*RSDSQ/(2*WSM) ENDIF IF(N1.GE.1) THEN !!! MAT(1,1)=MAT(1,1)+1/APR(1)**2 MAT(2,2)=MAT(2,2)+1/APR(2)**2 MAT(3,3)=MAT(3,3)+1/APR(3)**2 CALL INVERTN(3,MAT,IMAT) T1=T1+SNGL(IMAT(1,1)*W(1)+IMAT(1,2)*W(2)+IMAT(1,3)*W(3)) T2=T2+SNGL(IMAT(2,1)*W(1)+IMAT(2,2)*W(2)+IMAT(2,3)*W(3)) T3=T3+SNGL(IMAT(3,1)*W(1)+IMAT(3,2)*W(2)+IMAT(3,3)*W(3)) IF(MXSLP.NE.0) THEN Z2=T1**2+T2**2 IF(Z2.GT.MXSLP**2) THEN Z2=MXSLP/SQRT(Z2) T1=SNGL(Z2)*T1 T2=SNGL(Z2)*T2 ENDIF ENDIF TFND=.TRUE. RETURN ENDIF T1=0 T2=0 T3=0 TFND=.FALSE. RSDSQ=0 WSM=0 RETURN END