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$ Abstract
C
C 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 relative
C albedo, we are minimizing:
C
C ∑W(k)(O - P - dT1∂P/∂T1 - dT2∂P/∂T2 - dT3∂P/∂T3)^2
C k
C
C where the sum is over images, O is the observed brightness at a pixel and
C P is the predicted brightness. The weights W(k) and residuals O-P are
C carried in from the FIND_SLOPES subroutine. The partials with respect to
C T1 and T2 are determined numerically.
C
C There is no reason why these two subroutines could not be combined, but
C for the moment we are taking an "If it ain't broke ..." attitude.
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
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/O
C
C Variable I/O Description
C -------- --- --------------------------------------------------
C PMX
C NPIX
C RSD
C SP
C CP
C APR
C LAMBDA
C WGT
C T1
C T2
C T3
C USE
C PUSE
C TFND
C RSDSQ
C WSM
C S2S
C
C$ File_I/O
C
C Filename I/O Description
C ---------------------------- --- -------------------------------
C INIT_LITHOS.TXT I Setup file for the SPC toolkit
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 None
C
C$ SPC_subroutines_called
C INVERTN
C
C$ SPICELIB_functions_called
C None
C
C$ SPICELIB_subroutines_called
C None
C
C$ Called_by_SPC_Programs
C LITHOS
C SUBROUTINE FIND_SLOPES
C
% 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