C  PROCEDURE MOSAICX
C gfortran rawMOSAICXjrw.f /usr/local/toolkit/lib/spicelib.a
C 01 Nov 2018 - Added scaling, write out as TXT instead of pgm
C 05 Nov 2018 - EXTRACT_DATA no longer scales the values to go from 0 to 1.


      IMPLICIT NONE

      INTEGER               NTMP 
      PARAMETER            (NTMP=2500)        

      DOUBLE PRECISION      V(3)
      DOUBLE PRECISION      UX(3)
      DOUBLE PRECISION      UY(3)
      DOUBLE PRECISION      UZ(3)
      DOUBLE PRECISION      SCALE
      DOUBLE PRECISION      DNRAW(5001,5001)
      
      REAL*4                DN1(5001,5001)
      REAL*4                DNX(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                HT(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                ALB(-NTMP:NTMP,-NTMP:NTMP)
 
      INTEGER               I
      INTEGER               J
      INTEGER               QSZ
      INTEGER               NPX
      INTEGER               NLN

      INTEGER*2             dn(500000000)

      real                  scalingFactor
      real                  scalingTerm


      CHARACTER*6           MAPNM
      CHARACTER*12          PICNM
      CHARACTER*72          LMRKFILE
      CHARACTER*72          INFILE
      CHARACTER*72          OUTFILE
      CHARACTER*5001        CLINE
      
      LOGICAL               HUSE(-NTMP:NTMP,-NTMP:NTMP)
      LOGICAL               DONE(-NTMP:NTMP,-NTMP:NTMP)
      
      QSZ=NTMP
      DO I=-QSZ,QSZ
      DO J=-QSZ,QSZ
        HT(I,J)=0
        HUSE(I,J)=.FALSE.
        DN1(I+QSZ+1,J+QSZ+1)=0
        DNRAW(I+QSZ+1,J+QSZ+1)=0
        DONE(I,J)=.FALSE.
      ENDDO
      ENDDO
 
      WRITE(6,*) 'input 6-character map name.'
      READ(5,FMT='(A6)') MAPNM
      WRITE(6,*) 'Scaling Factor (multiply)'
      READ(5,*) scalingFactor
      WRITE(6,*) 'Offset (constant to add, 0 for if other for phase)'
      READ(5,*) scalingTerm

      LMRKFILE='MAPFILES/'//MAPNM//'.MAP'
      CALL READ_MAP(LMRKFILE,NTMP,QSZ,SCALE,V,UX,UY,UZ,HT,ALB)
      DO I=-QSZ,QSZ
      DO J=-QSZ,QSZ
        IF(ALB(I,J).GT.(0.005)) HUSE(I,J)=.TRUE.
      ENDDO
      ENDDO

      OPEN(UNIT=40,FILE='MOSAICX.TXT',STATUS='OLD')
10      CONTINUE
        READ(40,FMT='(A12)') PICNM
        IF(PICNM(1:3).NE.'END') THEN
          CALL PICINPT(PICNM,dn)
          CALL EXTRACT_DATAX(QSZ,SCALE,PICNM,
     .                        V,UX,UY,UZ,HUSE,HT,dn, dnx)
C Bob is changing the mapping from -qsz,qsz to 2*qsz+1 for DN1.
          DO I=-QSZ,QSZ
          DO J=-QSZ,QSZ
          IF((.NOT.DONE(I,J)).AND.(DNX(I,J).GT.0)) THEN
            DN1(I+QSZ+1,J+QSZ+1)=DNX(I,J)
            DONE(I,J)=.TRUE.
          ENDIF
          ENDDO
          ENDDO
          GO TO 10
        ENDIF
      CLOSE(UNIT=40)

      NPX=2*QSZ+1
      NLN=2*QSZ+1
      open(unit=10, file='TEMPFILE.gray', access='direct', 
     .     recl=NPX, status='UNKNOWN')
        DO J=1,NLN
          DO I=1,NPX
            CLINE(I:I)=CHAR(NINT(255*DN1(I,J)))
          ENDDO
          WRITE(10,REC=J) CLINE(1:NPX)
        ENDDO
      CLOSE(UNIT=10)
      infile= 'TEMPFILE.gray' 
      outfile= 'MOSAIC.pgm'  
      call raw2pgm(infile,outfile,npx,nln)
      open(unit=63, file=infile, status='old')
      close(unit=63, status='delete')
      WRITE(6,*)
      WRITE(6,*) 'gc MOSAIC.pgm'

C Write out raw values

      if (scalingTerm .eq. 0) then
        write (*,*) "# I/F scaling, no extra term"
        DO J=1,NLN
          DO I=1,NPX
            DNRAW(I,J)=DN1(I,J)*scalingFactor + scalingTerm
          ENDDO
        ENDDO
      else
        write (*,*) "# Phase scaling, with a x2"
        DO J=1,NLN
          DO I=1,NPX
            DNRAW(I,J)=2*DN1(I,J)*scalingFactor + scalingTerm
          ENDDO
        ENDDO
      endif

      infile= 'raw.txt'
      open(unit=12, file=infile, status='UNKNOWN')
        do j=1,2*QSZ+1
          do i=1,2*QSZ+1
           write (12, "(f12.5)", advance='no'),DNRAW(i, j)
          enddo
          write (12,*)
        enddo
      close(unit=12)


      STOP
      END

C$Procedure
 
      SUBROUTINE FIND_PICFILE(PICNM,PICFILE,EX)

C$ Abstract
C     This subroutine will search through directories
C         IMAGEFILES, IMAGEFILES1, IMAGEFILES2, IMAGEFILES3,  
C         IMAGEFILES4, IMAGEFILES5, IMAGEFILES6, IMAGEFILES7,
C         IMAGEFILES8, IMAGEFILES9 for the image file <PICNM> 
C         and returns the full path name to the file <PICFILE> and 
C         the logical variable <EX> indicating if it exists or not.
C
C$ Disclaimer
C     None
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*4             I, J
      INTEGER*4             SLEN
      CHARACTER*12          PICNM
      CHARACTER*72          PICFILE
      CHARACTER*72          ALTDIR
      CHARACTER*80          LINE
      LOGICAL               EX

      INTEGER*4             IFF
      SAVE                  IFF, ALTDIR
      DATA                  IFF/0/
 
C$ Variable_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     PICNM      I   Name of picture image file to search for
C     PICFILE    O   Path and file name link to the picture image file
C                    if files exists
C     EX         O   Logical variable indicating if file exists or not
C
C$ File_I/O
C
C     Filename                      I/O  Description
C     ----------------------------  ---  -------------------------------
C     None  
C
C$ Restrictions
C     None
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     SLEN
C
C$SPC_subroutines_called
C     None
C
C$SPICELIB_functions_called
C     None
C
C$SPICELIB_subroutines_called
C     None
C
C$ Called_by_SPC_Programs
C     AUTOREGISTER
C     AUTOREGISTERP
C     BLEMISHES
C     LIMBER
C     LITHOS
C     LITHOSP
C     RESIDUALS
C     SPHEREMAPSA
C     SUBROUTINE EXTRACT_DATA_PIC
C     SUBROUTINE FIND_PICFILE
C     SUBROUTINE LIMB_HEIGHTS
C     SUBROUTINE PICINPT
C
C$ History
C
C     2018_02_27:  A new INIT_LITHOS option ALTDIR= to allow IMAGEFILES 
C                  directories to be in other than the working directory 
C                  eg ALTDIR='../ROSETTA'
C 
C

      IF(IFF.EQ.0) THEN
        IFF=1
        ALTDIR='NONE'
        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:7).EQ.'ALTDIR=') READ(LINE(8:80),*) ALTDIR 
          GO TO 13 
        ENDIF
      ENDIF
      CLOSE(UNIT=25)

      I=SLEN(PICNM)
      PICFILE='./IMAGEFILES/'//PICNM(1:I)//'.DAT'
      INQUIRE(FILE=PICFILE,EXIST=EX)
      IF(EX) RETURN
      PICFILE='./IMAGEFILES1/'//PICNM(1:I)//'.DAT'
      INQUIRE(FILE=PICFILE,EXIST=EX)
      IF(EX) RETURN
      PICFILE='./IMAGEFILES2/'//PICNM(1:I)//'.DAT'
      INQUIRE(FILE=PICFILE,EXIST=EX)
      IF(EX) RETURN
      PICFILE='./IMAGEFILES3/'//PICNM(1:I)//'.DAT'
      INQUIRE(FILE=PICFILE,EXIST=EX)
      IF(EX) RETURN
      PICFILE='./IMAGEFILES4/'//PICNM(1:I)//'.DAT'
      INQUIRE(FILE=PICFILE,EXIST=EX)
      IF(EX) RETURN
      PICFILE='./IMAGEFILES5/'//PICNM(1:I)//'.DAT'
      INQUIRE(FILE=PICFILE,EXIST=EX)
      IF(EX) RETURN
      PICFILE='./IMAGEFILES6/'//PICNM(1:I)//'.DAT'
      INQUIRE(FILE=PICFILE,EXIST=EX)
      IF(EX) RETURN
      PICFILE='./IMAGEFILES7/'//PICNM(1:I)//'.DAT'
      INQUIRE(FILE=PICFILE,EXIST=EX)
      IF(EX) RETURN
      PICFILE='./IMAGEFILES8/'//PICNM(1:I)//'.DAT'
      INQUIRE(FILE=PICFILE,EXIST=EX)
      IF(EX) RETURN
      PICFILE='./IMAGEFILES9/'//PICNM(1:I)//'.DAT'
      INQUIRE(FILE=PICFILE,EXIST=EX)
      IF(EX) RETURN
      IF(ALTDIR(1:4).NE.'NONE') THEN
        J=SLEN(ALTDIR)
        PICFILE=ALTDIR(1:J)//'/IMAGEFILES/'//PICNM(1:I)//'.DAT'
        INQUIRE(FILE=PICFILE,EXIST=EX)
        IF(EX) RETURN
        PICFILE=ALTDIR(1:J)//'/IMAGEFILES1/'//PICNM(1:I)//'.DAT'
        INQUIRE(FILE=PICFILE,EXIST=EX)
        IF(EX) RETURN
        PICFILE=ALTDIR(1:J)//'/IMAGEFILES2/'//PICNM(1:I)//'.DAT'
        INQUIRE(FILE=PICFILE,EXIST=EX)
        IF(EX) RETURN
        PICFILE=ALTDIR(1:J)//'/IMAGEFILES3/'//PICNM(1:I)//'.DAT'
        INQUIRE(FILE=PICFILE,EXIST=EX)
        IF(EX) RETURN
        PICFILE=ALTDIR(1:J)//'/IMAGEFILES4/'//PICNM(1:I)//'.DAT'
        INQUIRE(FILE=PICFILE,EXIST=EX)
        IF(EX) RETURN
        PICFILE=ALTDIR(1:J)//'/IMAGEFILES5/'//PICNM(1:I)//'.DAT'
        INQUIRE(FILE=PICFILE,EXIST=EX)
        IF(EX) RETURN
        PICFILE=ALTDIR(1:J)//'/IMAGEFILES6/'//PICNM(1:I)//'.DAT'
        INQUIRE(FILE=PICFILE,EXIST=EX)
        IF(EX) RETURN
        PICFILE=ALTDIR(1:J)//'/IMAGEFILES7/'//PICNM(1:I)//'.DAT'
        INQUIRE(FILE=PICFILE,EXIST=EX)
        IF(EX) RETURN
        PICFILE=ALTDIR(1:J)//'/IMAGEFILES8/'//PICNM(1:I)//'.DAT'
        INQUIRE(FILE=PICFILE,EXIST=EX)
        IF(EX) RETURN
        PICFILE=ALTDIR(1:J)//'/IMAGEFILES9/'//PICNM(1:I)//'.DAT'
        INQUIRE(FILE=PICFILE,EXIST=EX)
        IF(EX) RETURN
      ENDIF
      PICFILE='NONE'

      RETURN
      END

c  ......................................................
      function brt0(z,npx,nln,dn)
c  ......................................................
      
      IMPLICIT NONE

      real*8        z(2), x, y, brt0
      integer*4     i, j, b0, b1, b2, b3, npx, nln
      integer*2     dn(500000000)

      brt0=0
      if((z(1).ge.1).and.(z(1).lt.npx).and.
     .   (z(2).ge.1).and.(z(2).lt.nln)) then
        i=int(z(1))
        j=int(z(2))
        x=z(1)-i
        y=z(2)-j
        b0=dn(i+npx*(j-1))
        b1=dn(i+1+npx*(j-1))
        b2=dn(i+npx*j)
        b3=dn(i+1+npx*j)
        if((b0.eq.0).and.(b1*b2*b3.ne.0)) b0=max(1,b1+b2-b3)
        if((b1.eq.0).and.(b2*b3*b0.ne.0)) b1=max(1,b0+b3-b2)
        if((b2.eq.0).and.(b3*b0*b1.ne.0)) b2=max(1,b3+b0-b1)
        if((b3.eq.0).and.(b0*b1*b2.ne.0)) b3=max(1,b2+b1-b0)
        if(b0*b1*b2*b3.eq.0) return
        b3=b0-b1-b2+b3
        b1=b1-b0
        b2=b2-b0
        brt0=b0+b1*x+b2*y+b3*x*y
        brt0=max(0.d0,brt0)
        return
      endif

      return
      end

c   ................................................
      SUBROUTINE V2IMGPL(V,V0,PICNM,NPX,NLN,MMFL,CTR,KMAT,D,
     .                   CX,CY,CZ, USE,IMGPL)
c   ................................................

      IMPLICIT NONE

      DOUBLE PRECISION      VDOT
      DOUBLE PRECISION      VNORM

      DOUBLE PRECISION      IMGPL(2)
      DOUBLE PRECISION      MMFL
      DOUBLE PRECISION      KMAT(2,3)
      DOUBLE PRECISION      D(4)
      DOUBLE PRECISION      CTR(2)
      DOUBLE PRECISION      MM(2)
      DOUBLE PRECISION      Z(2)
      DOUBLE PRECISION      V0(3)
      DOUBLE PRECISION      V(3)
      DOUBLE PRECISION      V1(3)
      DOUBLE PRECISION      U(3)
      DOUBLE PRECISION      CX(3)      
      DOUBLE PRECISION      CY(3)
      DOUBLE PRECISION      CZ(3)
      DOUBLE PRECISION      WX(3)
      DOUBLE PRECISION      WY(3)
      DOUBLE PRECISION      WZ(3)
      DOUBLE PRECISION      Z1
      DOUBLE PRECISION      Z2
      DOUBLE PRECISION      Z3
      DOUBLE PRECISION      R0(3)
      DOUBLE PRECISION      U0(3)
      DOUBLE PRECISION      CAX(20,3)
      DOUBLE PRECISION      TRJ_AXIS(3)
      DOUBLE PRECISION      TRJ(3)
      DOUBLE PRECISION      POL_AXIS(3)
      DOUBLE PRECISION      POL
      DOUBLE PRECISION      X0, X1, Y0, Y1
          
      INTEGER               I
      INTEGER               K
      INTEGER               KMAX
      INTEGER               NPX
      INTEGER               NLN
      INTEGER               SLEN
          
      CHARACTER*12          PICNM
      CHARACTER*72          PICTFILE
      CHARACTER*80          LINE
          
      LOGICAL               USE
      LOGICAL               LS

      USE=.FALSE.
      LS=.FALSE.
      I=SLEN(PICNM)
      PICTFILE='SUMFILES/'//PICNM(1:I)//'.SUM'
      OPEN(UNIT=32,FILE=PICTFILE,STATUS='OLD')
        READ(32,FMT='(A80)') LINE
        IF(LINE(14:14).EQ.'L') THEN
          LS=.TRUE.
10        CONTINUE
          READ(32,FMT='(A80)') LINE
          IF(LINE(1:14).NE.'LINE SCAN DATA') GO TO 10
          READ(32,*) (U0(I), I=1,3) 
          READ(32,*) (R0(I), I=1,3)
          READ(32,*) (POL_AXIS(I), I=1,3)
          READ(32,*) POL
          READ(32,*) (TRJ_AXIS(I), I=1,3)
          READ(32,*) (TRJ(I), I=1,3)
          READ(32,*) KMAX
          DO K=1,KMAX
            READ(32,*) (CAX(K,I), I=1,3)
          ENDDO
        ENDIF
      CLOSE(UNIT=32)

      IF(LS) THEN
        R0(1)=VNORM(V0)
        CALL VHAT(V0,U0)
        Z(2)=CTR(2)
        DO K=1,3
          X0=Z(2)-50
          CALL SCOBJ(R0,U0,TRJ_AXIS,TRJ,POL_AXIS,POL,CTR(2),X0,V1)
          CALL POINTING(CX,CY,CZ,CAX,KMAX,TRJ_AXIS,TRJ,
     .                      POL_AXIS,POL,CTR(2),X0,WX,WY,WZ)
          CALL VADD(V1,V,U)
          Z1=VDOT(U,WX)
          Z2=VDOT(U,WY)
          Z3=VDOT(U,WZ)
          call getmm(picnm,mmfl,z1,z2,z3,mm)
          Y0=MM(2)
          X1=Z(2)+50
          CALL SCOBJ(R0,U0,TRJ_AXIS,TRJ,POL_AXIS,POL,CTR(2),X1,V1)
          CALL POINTING(CX,CY,CZ,CAX,KMAX,TRJ_AXIS,TRJ,
     .                      POL_AXIS,POL,CTR(2),X1,WX,WY,WZ)
          CALL VADD(V1,V,U)
          Z1=VDOT(U,WX)
          Z2=VDOT(U,WY)
          Z3=VDOT(U,WZ)
          call getmm(picnm,mmfl,z1,z2,z3,mm)
          Y1=MM(2)
          Z(2)=(X0*Y1-X1*Y0)/(Y1-Y0)
        ENDDO
        X0=Z(2)
        CALL SCOBJ(R0,U0,TRJ_AXIS,TRJ,POL_AXIS,POL,CTR(2),X0,V1)
        CALL POINTING(CX,CY,CZ,CAX,KMAX,TRJ_AXIS,TRJ,
     .                    POL_AXIS,POL,CTR(2),X0,WX,WY,WZ)
        CALL VADD(V1,V,U)
        Z1=VDOT(U,WX)
        Z2=VDOT(U,WY)
        Z3=VDOT(U,WZ)
        call getmm(picnm,mmfl,z1,z2,z3,mm)
        IMGPL(1)=CTR(1)+KMAT(1,1)*MM(1)
        IMGPL(2)=X0
        IF((IMGPL(1).GE.1).AND.(IMGPL(1).LE.NPX).AND.
     .     (IMGPL(2).GE.1).AND.(IMGPL(2).LE.NLN)) THEN
          USE=.TRUE.
        ENDIF
        RETURN
      ENDIF

      CALL VADD(V0,V,U)
      CALL VHAT(U,U)
      Z1=VDOT(U,CX)
      Z2=VDOT(U,CY)
      Z3=VDOT(U,CZ)
      call getmm(picnm,mmfl,z1,z2,z3,mm)
      CALL MMPX(KMAT,CTR,D,MM,IMGPL)
      IF((IMGPL(1).GE.1).AND.(IMGPL(1).LE.NPX).AND.
     .   (IMGPL(2).GE.1).AND.(IMGPL(2).LE.NLN)) THEN
        USE=.TRUE.
      ENDIF

      RETURN
      END

c   ................................................
      SUBROUTINE SCOBJ(R0,U0,TRJ_AXIS,TRJ,POL_AXIS,POL,LN0,LN,V)
c   ................................................

      IMPLICIT NONE

      DOUBLE PRECISION  R0(3)
      DOUBLE PRECISION  U0(3)
      DOUBLE PRECISION  TRJ_AXIS(3)
      DOUBLE PRECISION  TRJ(3)
      DOUBLE PRECISION  POL_AXIS(3)
      DOUBLE PRECISION  POL
      DOUBLE PRECISION  LN0
      DOUBLE PRECISION  LN
      DOUBLE PRECISION  V(3)
      DOUBLE PRECISION  Z

      Z=TRJ(1)*(LN-LN0)+TRJ(2)*(LN-LN0)**2+TRJ(3)*(LN-LN0)**3
      CALL VROTV(U0,TRJ_AXIS,Z,V)
      Z=R0(1)+R0(2)*(LN-LN0)+R0(3)*(LN-LN0)**2
      CALL VSCL(Z,V,V)
      Z=-POL*(LN-LN0)
      CALL VROTV(V,POL_AXIS,Z,V)
          
      RETURN
      END

c   ................................................
      SUBROUTINE POINTING(CX,CY,CZ,CAX,IMAX,TRJ_AXIS,TRJ,
     .                    POL_AXIS,POL,LN0,LN,WX,WY,WZ)
c   ................................................

      IMPLICIT NONE

      DOUBLE PRECISION  CX(3)
      DOUBLE PRECISION  CY(3)
      DOUBLE PRECISION  CZ(3)
      DOUBLE PRECISION  CAX(20,3)
      DOUBLE PRECISION  TRJ_AXIS(3)
      DOUBLE PRECISION  TRJ(3)
      DOUBLE PRECISION  POL_AXIS(3)
      DOUBLE PRECISION  POL
      DOUBLE PRECISION  LN0
      DOUBLE PRECISION  LN
      DOUBLE PRECISION  WX(3)
      DOUBLE PRECISION  WY(3)
      DOUBLE PRECISION  WZ(3)
      DOUBLE PRECISION  V(3)
      DOUBLE PRECISION  CMAT(3,3)
      DOUBLE PRECISION  XMAT(3,3)
      DOUBLE PRECISION  Z
      
      INTEGER           I
      INTEGER           J
      INTEGER           IMAX

      DO J=1,3
        CMAT(1,J)=CX(J)
        CMAT(2,J)=CY(J)
        CMAT(3,J)=CZ(J)
      ENDDO
      IF(LN.NE.LN0) THEN
        DO J=1,3
          V(J)=0
          DO I=1,IMAX
            V(J)=V(J)+CAX(I,J)*(LN-LN0)**I
          ENDDO
        ENDDO
        CALL ROT2MAT(V,XMAT)
        CALL MXM(XMAT,CMAT,CMAT)
      ENDIF
      DO J=1,3
        WX(J)=CMAT(1,J)
        WY(J)=CMAT(2,J)
        WZ(J)=CMAT(3,J)
      ENDDO
      Z=TRJ(1)*(LN-LN0)+TRJ(2)*(LN-LN0)**2+TRJ(3)*(LN-LN0)**3
      CALL VROTV(WX,TRJ_AXIS,Z,WX)
      CALL VROTV(WY,TRJ_AXIS,Z,WY)
      CALL VROTV(WZ,TRJ_AXIS,Z,WZ)
      Z=-POL*(LN-LN0)
      CALL VROTV(WX,POL_AXIS,Z,WX)
      CALL VROTV(WY,POL_AXIS,Z,WY)
      CALL VROTV(WZ,POL_AXIS,Z,WZ)
          
      RETURN
      END

c   ................................................
      SUBROUTINE ROT2MAT(V,M)
c   ................................................

      implicit none

      real*8     V(3), M(3,3), U(3), Z, VNORM
      integer*4  i, j
          
      Z=VNORM(V)
      IF(Z.EQ.0) THEN
        DO I=1,3
        DO J=1,3
          M(I,J)=0.D0
        ENDDO
        ENDDO
        DO I=1,3
          M(I,I)=1.D0
        ENDDO
        RETURN
      ENDIF
      CALL VHAT(V,U)
            
      do i=1,3
      do j=1,3
        M(i,j)=u(i)*u(j)*(1.d0-cos(z))
      enddo
      enddo
      M(1,1)=M(1,1)+cos(z)
      M(2,2)=M(2,2)+cos(z)
      M(3,3)=M(3,3)+cos(z)
      M(1,2)=M(1,2)-u(3)*sin(z)
      M(2,1)=M(2,1)+u(3)*sin(z)
      M(2,3)=M(2,3)-u(1)*sin(z)
      M(3,2)=M(3,2)+u(1)*sin(z)
      M(3,1)=M(3,1)-u(2)*sin(z)
      M(1,3)=M(1,3)+u(2)*sin(z)
        
      return
      end

c   ................................................
      SUBROUTINE MAT2ROT(M,V)
c   ................................................

      IMPLICIT NONE

      DOUBLE PRECISION  M(3,3)
      DOUBLE PRECISION  V(3)
      DOUBLE PRECISION  U(3)
      DOUBLE PRECISION  VDOT
      DOUBLE PRECISION  PI
      DOUBLE PRECISION  Z
      DOUBLE PRECISION  X
      DOUBLE PRECISION  W(3)

      INTEGER           I
      INTEGER           J

      W(1)=(M(2,2)-1.D0)*(M(3,3)-1.D0)-M(2,3)*M(3,2)
      W(2)=(M(3,3)-1.D0)*(M(1,1)-1.D0)-M(3,1)*M(1,3)
      W(3)=(M(1,1)-1.D0)*(M(2,2)-1.D0)-M(1,2)*M(2,1)

      X=0
      DO I=1,3
        IF(ABS(W(I)).GT.X) THEN
          X=ABS(W(I))
          J=I
        ENDIF
      ENDDO

      IF(J.EQ.1) THEN
        U(1)=1.D0
        U(2)=(M(2,3)*M(3,1)-M(2,1)*(M(3,3)-1.D0))/W(1)
        U(3)=(M(3,2)*M(2,1)-M(3,1)*(M(2,2)-1.D0))/W(1)
      ENDIF
      IF(J.EQ.2) THEN
        U(2)=1.D0
        U(3)=(M(3,1)*M(1,2)-M(3,2)*(M(1,1)-1.D0))/W(2)
        U(1)=(M(1,3)*M(3,2)-M(1,2)*(M(3,3)-1.D0))/W(2)
      ENDIF
      IF(J.EQ.3) THEN
        U(3)=1.D0
        U(1)=(M(1,2)*M(2,3)-M(1,3)*(M(2,2)-1.D0))/W(3)
        U(2)=(M(2,1)*M(1,3)-M(2,3)*(M(1,1)-1.D0))/W(3)
      ENDIF
      CALL VHAT(U,U)

      X=2
      DO I=1,3
        W(I)=0
        IF(ABS(U(I)).LT.X) THEN
          X=ABS(U(I))
          J=I
        ENDIF
      ENDDO
      W(J)=1
      CALL MXV(M,W,V)
      X=(VDOT(W,V)-VDOT(W,U)**2)/(1.D0-VDOT(W,U)**2)
      CALL VCRSS(U,V,V)
      Z=ASIN(VDOT(W,V)/(VDOT(W,U)**2-1.D0))
      IF(X.LT.0) THEN
        Z=PI()-Z
      ENDIF
      CALL VSCL(Z,U,V)
          
      RETURN  
      END
          
C   ..................................................
      SUBROUTINE MMPX(KMAT,CTR,D,MM,PX)
C   ..................................................

      IMPLICIT NONE
          
      DOUBLE PRECISION      KMAT(2,3)
      DOUBLE PRECISION      CTR(2)
      DOUBLE PRECISION      MM(2)  
      DOUBLE PRECISION      PX(2)
      DOUBLE PRECISION      D(4)
      DOUBLE PRECISION      X
      DOUBLE PRECISION      Y
      DOUBLE PRECISION      Z
      DOUBLE PRECISION      RSQ
      INTEGER               K
      LOGICAL               USE
      
      X=MM(1)*(1.d0+D(1))
      Y=MM(2)*(1.d0+D(1))

      USE=.FALSE.
      DO K=2,4
        USE=USE.OR.(D(K).NE.0)
      ENDDO
      
      IF(USE) THEN
        RSQ=X**2+Y**2
        Z=1.d0+D(2)*RSQ+D(3)*Y+D(4)*X
        X=X*Z
        Y=Y*Z
      ENDIF
          
      PX(1)=CTR(1)+KMAT(1,1)*X+KMAT(1,2)*Y+KMAT(1,3)*X*Y
      PX(2)=CTR(2)+KMAT(2,1)*X+KMAT(2,2)*Y+KMAT(2,3)*X*Y
          
      RETURN
      END

C   ..................................................
      subroutine getmm(picnm,mmfl,z1,z2,z3,z)
C   ..................................................

      IMPLICIT NONE

      double precision      mmfl
      double precision      z1, z2, z3
      double precision      z(2)
      double precision      r
      double precision      pd(0:100,100)
      double precision      X(100)

      INTEGER               JD
      INTEGER               KD
      INTEGER               I
      INTEGER               J

      CHARACTER*12          PICNM
      CHARACTER*12          TMPL(100)
      CHARACTER*72          TYPE(100)
      CHARACTER*80          LINE

      LOGICAL               EX

      INTEGER*4             IFF
      SAVE                  IFF, KD, TMPL, TYPE, PD
      data                  iff/0/

      IF(IFF.EQ.0) THEN
        IFF=1
        KD=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:8).EQ.'DISTORT=') THEN
                KD=KD+1
                READ(LINE(9:80),*) TMPL(KD), TYPE(KD), J
                READ(25,*) (PD(I,KD), I=0,J)
              ENDIF
              GO TO 13
            ENDIF 
          CLOSE(UNIT=25)
        ENDIF
      ENDIF

      IF(KD.NE.0) THEN
        JD=0
        DO J=1,KD
          DO I=1,12
            IF((TMPL(J)(I:I).NE.'#').AND.(TMPL(J)(I:I).NE.' ')
     .         .AND.(PICNM(I:I).NE.TMPL(J)(I:I))) GO TO 10
          ENDDO
          JD=J
10        CONTINUE
        ENDDO
        IF(JD.NE.0) THEN
          IF(TYPE(JD).EQ.'CUBIC') THEN
            z(1)=(1+PD(0,JD))*mmfl*z1/z3
            z(2)=(1+PD(0,JD))*mmfl*z2/z3
            X( 1)=1.D0
            X( 2)=z(1)
            X( 3)=z(2)            
            X( 4)=z(1)*z(1)
            X (5)=z(1)*z(2)            
            X( 6)=z(2)*z(2)            
            X( 7)=z(1)*z(1)*z(1)
            X( 8)=z(1)*z(1)*z(2)            
            X( 9)=z(1)*z(2)*z(2)            
            X(10)=z(2)*z(2)*z(2)   
            z(1)=0
            z(2)=0  
            DO I=1,10
              z(1)=z(1)+PD(I,JD)*X(I)
              z(2)=z(2)+PD(I+10,JD)*X(I)
            ENDDO
            RETURN
          ENDIF         
          IF(TYPE(JD).EQ.'RADIAL') THEN
            z(1)=(1+PD(0,JD))*mmfl*z1/z3
            z(2)=(1+PD(0,JD))*mmfl*z2/z3
            X(1)=z(1)**2+z(2)**2
            z(1)=z(1)*(1+PD(1,JD)*X(1)+PD(2,JD)*X(1)**2)           
            z(2)=z(2)*(1+PD(1,JD)*X(1)+PD(2,JD)*X(1)**2)           
            RETURN
          ENDIF         
          IF(TYPE(JD).EQ.'OWEN') THEN
            z(1)=(1+PD(0,JD))*mmfl*z1/z3
            z(2)=(1+PD(0,JD))*mmfl*z2/z3
            R=SQRT(z(1)**2+z(2)**2)
            X(1)=z(1)
            X(2)=z(2) 
            z(1)=z(1)-PD(1,JD)*X(2)*R+PD(2,JD)*X(1)*R**2
     .               -PD(3,JD)*X(2)*R**3+PD(4,JD)*X(1)*R**4
     .               +PD(5,JD)*X(1)*X(2)+PD(6,JD)*X(1)*X(1)
            z(2)=z(2)+PD(1,JD)*X(1)*R+PD(2,JD)*X(2)*R**2
     .               +PD(3,JD)*X(1)*R**3+PD(4,JD)*X(2)*R**4
     .               +PD(5,JD)*X(2)*X(2)+PD(6,JD)*X(1)*X(2)
            RETURN
          ENDIF         
          IF(TYPE(JD).EQ.'LINESCAN') THEN
            z(2)=mmfl*z2/z3
            Z(1)=(1+PD(0,JD))*mmfl*z1/z3
            R=SQRT(Z(1)**2+PD(1,JD)**2)
            Z(1)=Z(1)*(1+PD(2,JD)*R**2+PD(3,JD)*R**3)
            RETURN
          ENDIF         
          IF(TYPE(JD).EQ.'APOLLO') THEN
            z(1)=(1+PD(0,JD))*mmfl*z1/z3
            z(2)=(1+PD(0,JD))*mmfl*z2/z3
            X(1)=z(1)**2+z(2)**2
            X(2)=X(1)*X(1)
            X(3)=X(1)*X(2)
            z(1)=z(1)*(1+PD(1,JD)*X(1)+PD(2,JD)*X(2)+PD(3,JD)*X(3))           
            z(2)=z(2)*(1+PD(1,JD)*X(1)+PD(2,JD)*X(2)+PD(3,JD)*X(3))
            X(4)=PD(4,JD)*X(1)+PD(5,JD)*X(2)
            Z(1)=Z(1)+X(4)*SIN(PD(6,JD))           
            Z(2)=Z(2)-X(4)*COS(PD(6,JD))           
            RETURN
          ENDIF         
          IF(TYPE(JD).EQ.'HIRISE') THEN
            z(1)=(1+PD(0,JD))*mmfl*z1/z3
            z(2)=(1+PD(0,JD))*mmfl*z2/z3
            z(1)=z(1)-PD(1,JD)
            z(2)=z(2)-PD(2,JD)
            RETURN
          ENDIF         
        ENDIF
      ENDIF

      z(1)=mmfl*z1/z3
      z(2)=mmfl*z2/z3

      RETURN
      END
      
c  ...............................................................
      SUBROUTINE IMGPL2VN(PICNM,IMGPL,USE,V,N)
c  ...............................................................

      IMPLICIT NONE

      DOUBLE PRECISION      V0(3)
      DOUBLE PRECISION      W0(3)
      DOUBLE PRECISION      V(3)
      DOUBLE PRECISION      CX(3)      
      DOUBLE PRECISION      CY(3)
      DOUBLE PRECISION      CZ(3)
      DOUBLE PRECISION      WX(3)      
      DOUBLE PRECISION      WY(3)
      DOUBLE PRECISION      WZ(3)
      DOUBLE PRECISION      N(3)
      DOUBLE PRECISION      IMGPL(2)
      DOUBLE PRECISION      Z1, Z2, Z3
      DOUBLE PRECISION      MMFL
      DOUBLE PRECISION      KMAT(2,3)
      DOUBLE PRECISION      D(4)
      DOUBLE PRECISION      CTR(2)
      DOUBLE PRECISION      MM(2)
      DOUBLE PRECISION      R0(3)
      DOUBLE PRECISION      U0(3)
      DOUBLE PRECISION      CAX(20,3)
      DOUBLE PRECISION      TRJ_AXIS(3)
      DOUBLE PRECISION      TRJ(3)
      DOUBLE PRECISION      POL_AXIS(3)
      DOUBLE PRECISION      POL

      INTEGER               NPX
      INTEGER               NLN
      INTEGER               KMAX
      INTEGER               SLEN
      INTEGER               I
      INTEGER               K

      CHARACTER*12          PICNM
      CHARACTER*72          PICTFILE
      character*80          LINE

      LOGICAL               USE
      LOGICAL               LS

      LS=.FALSE.
      D(1)=0.D0
      D(2)=0.D0
      D(3)=0.D0
      D(4)=0.D0
      I=SLEN(PICNM)
      PICTFILE='SUMFILES/'//PICNM(1:I)//'.SUM'
      OPEN(UNIT=10,FILE=PICTFILE,STATUS='OLD')
        READ(10,FMT='(A80)') LINE
        IF(LINE(14:14).EQ.'L')  LS=.TRUE.
        READ(10,*)
        READ(10,*) NPX, NLN
        READ(10,*) MMFL, CTR(1), CTR(2)
        READ(10,*) (V0(I), I=1,3) 
        READ(10,*) (CX(I), I=1,3)
        READ(10,*) (CY(I), I=1,3)
        READ(10,*) (CZ(I), I=1,3)
        READ(10,*)! SZ
        READ(10,*)  KMAT(1,1), KMAT(1,2), KMAT(1,3),
     .              KMAT(2,1), KMAT(2,2), KMAT(2,3)
        READ(10,FMT='(A80)') LINE
        IF(LINE(64:73).EQ.'DISTORTION') THEN
          READ(LINE,*) (D(I), I=1,4)
        ENDIF
        IF(LS) THEN
10        CONTINUE
          READ(10,FMT='(A80)') LINE
          IF(LINE(1:14).NE.'LINE SCAN DATA') GO TO 10
          READ(10,*) (U0(I), I=1,3) 
          READ(10,*) (R0(I), I=1,3)
          READ(10,*) (POL_AXIS(I), I=1,3)
          READ(10,*) POL
          READ(10,*) (TRJ_AXIS(I), I=1,3)
          READ(10,*) (TRJ(I), I=1,3)
          READ(10,*) KMAX
          DO K=1,KMAX
            READ(10,*) (CAX(K,I), I=1,3)
          ENDDO
        ENDIF
      CLOSE(UNIT=10)

      IF(LS) THEN
        CALL SCOBJ(R0,U0,TRJ_AXIS,TRJ,POL_AXIS,POL,CTR(2),IMGPL(2),W0)
        CALL POINTING(CX,CY,CZ,CAX,KMAX,TRJ_AXIS,TRJ,
     .                    POL_AXIS,POL,CTR(2),IMGPL(2),WX,WY,WZ)
        Z1=0
        Z2=0
        Z3=1
        call getmm(picnm,mmfl,z1,z2,z3,mm)
        MM(1)=-MM(1)+(IMGPL(1)-CTR(1))/KMAT(1,1)
        MM(2)=-MM(2)
        CALL MM2VN(MM,MMFL,W0,WX,WY,WZ,USE,V,N)
      ELSE
        CALL VEQU(V0,W0)
        CALL VEQU(CX,WX)
        CALL VEQU(CY,WY)
        CALL VEQU(CZ,WZ)
        CALL PXMM(PICNM,MMFL,KMAT,CTR,D,IMGPL,MM)
        CALL MM2VN(MM,MMFL,W0,WX,WY,WZ,USE,V,N)
      ENDIF

      return
      end

c  ...............................................................
      subroutine PICINPT(PICNM,dn)
c  ...............................................................
      
      IMPLICIT NONE

      DOUBLE PRECISION      Z1, Z2

      CHARACTER*3           ORDER
      CHARACTER*12          PICNM
      CHARACTER*72          PICTFILE
      CHARACTER*99999       INLINE
      
      INTEGER               NPX
      INTEGER               NLN
      INTEGER               I
      INTEGER               J
      INTEGER               I1
      INTEGER               J1
      INTEGER               I2
      INTEGER               J2
      INTEGER               SLEN
      INTEGER               BYTES
      INTEGER               d1(-1:1,-1:1)
      INTEGER               d2(-1:1,-1:1)

      INTEGER*2             dn(500000000)

      LOGICAL               EX

      I=SLEN(PICNM)      
      PICTFILE='SUMFILES/'//PICNM(1:I)//'.SUM'
      OPEN(UNIT=10,FILE=PICTFILE,STATUS='OLD')
        READ(10,*)
        READ(10,*)
        READ(10,*) NPX, NLN
      CLOSE(UNIT=10)
      CALL FIND_PICFILE(PICNM,PICTFILE,EX)
      IF(.NOT.EX) THEN
        DO J=1,NLN
        DO I=1,NPX
          DN(I+NPX*(J-1))=0
        ENDDO
        ENDDO
        RETURN
      ENDIF

      bytes=2
      open(unit=10, file=pictfile, access='direct',
     .     recl=2*npx, status='old')
        read(10,rec=nln, err=20) inline(1:2*npx)
        go to 21
20      bytes=1
21      continue
      close(unit=10)

      if(bytes.eq.2) then
        open(unit=10, file=pictfile, access='direct',
     .       recl=2*npx, status='old')
          z1=0
          z2=0
          do i1=1,99,2
          do j1=1,99,2
            i=(5+i1*npx)/100
            j=(5+j1*nln)/100
            do i2=-1,1
            do j2=-1,1
              read(10,rec=j+j2) inline(1:2*npx)
              d1(i2,j2)=ichar(inline(2*(i+i2)-1:2*(i+i2)-1)) 
              d2(i2,j2)=ichar(inline(2*(i+i2):2*(i+i2))) 
            enddo
            enddo
            z1=z1+(4*d1(0,0)-d1(-1,0)-d1(1,0)-d1(0,-1)-d1(0,1))**2
            z2=z2+(4*d2(0,0)-d2(-1,0)-d2(1,0)-d2(0,-1)-d2(0,1))**2
          enddo
          enddo 
        close(unit=10)
        if(z1.gt.z2) order='LSB'
        if(z2.gt.z1) order='MSB'
      endif

      if(bytes.eq.1) then
        open(unit=10, file=PICTFILE, access='direct', 
     .       recl=NPX, status='old')
          do j=1,NLN
            read(10,rec=j) inline(1:npx)
            do i=1,NPX
              DN(I+NPX*(J-1))=ichar(inline(I:I))
            enddo
          enddo
        close (unit=10)
      endif
      if(bytes.eq.2) then
        open(unit=10, file=PICTFILE, access='direct', 
     .       recl=2*NPX, status='old')
          do j=1,NLN
            read(10,rec=j) inline(1:2*npx)
            do i=1,NPX
              if(order.eq.'MSB') DN(I+NPX*(J-1))=
     .          ichar(inline(2*I:2*I))+256*ichar(inline(2*I-1:2*I-1))
              if(order.eq.'LSB') DN(I+NPX*(J-1))=
     .          256*ichar(inline(2*I:2*I))+ichar(inline(2*I-1:2*I-1))
            enddo
          enddo
        close(unit=10)
      endif

      return
      end

c  ...............................................................
      SUBROUTINE EXTRACT_DATAX(QSZ,SCALE,PICNM,
     .                        V,UX,UY,UZ,HUSE,HT,dn, dnx)
c  ...............................................................

      implicit none
      
      INTEGER               NTMP 
      PARAMETER            (NTMP=2500)        

      DOUBLE PRECISION      V(3)
      DOUBLE PRECISION      W(3)
      DOUBLE PRECISION      UX(3)
      DOUBLE PRECISION      UY(3)
      DOUBLE PRECISION      UZ(3)
      DOUBLE PRECISION      CX(3)
      DOUBLE PRECISION      CY(3)
      DOUBLE PRECISION      CZ(3)
      DOUBLE PRECISION      V0(3)
      DOUBLE PRECISION      SZ(3)
      DOUBLE PRECISION      MMFL
      DOUBLE PRECISION      KMAT(2,3)
      DOUBLE PRECISION      D(4)
      DOUBLE PRECISION      CTR(2)
      DOUBLE PRECISION      Z(2)
      DOUBLE PRECISION      Z2
      DOUBLE PRECISION      BRT0
      DOUBLE PRECISION      SCALE
      
      REAL*4                REAL4
      REAL*4                DNX(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                HT(-NTMP:NTMP,-NTMP:NTMP)
          
      INTEGER               I
      INTEGER               J
      INTEGER               QSZ
      INTEGER               NPX
      INTEGER               NLN
      INTEGER               T1
      INTEGER               T2
      INTEGER               SLEN

      INTEGER*2             dn(500000000)

      CHARACTER*12          PICNM         
      CHARACTER*72          PICTFILE
      CHARACTER*80          LINE

      LOGICAL               HUSE(-NTMP:NTMP,-NTMP:NTMP)
      LOGICAL               EX

      I=SLEN(PICNM)
      PICTFILE='SUMFILES/'//PICNM(1:I)//'.SUM'
      OPEN(UNIT=20,FILE=PICTFILE,STATUS='OLD')
        READ(20,*)
        READ(20,*)
        READ(20,*)  NPX, NLN, T1, T2
        READ(20,*)  MMFL,CTR(1), CTR(2)
        READ(20,*) (V0(I), I=1,3) 
        READ(20,*) (CX(I), I=1,3)
        READ(20,*) (CY(I), I=1,3)
        READ(20,*) (CZ(I), I=1,3)
        READ(20,*) (SZ(I), I=1,3)
        READ(20,*)  KMAT(1,1), KMAT(1,2), KMAT(1,3),
     .              KMAT(2,1), KMAT(2,2), KMAT(2,3)
        READ(20,FMT='(A80)') LINE
        IF(LINE(64:73).EQ.'DISTORTION') THEN
          READ(LINE,*) (D(I), I=1,4)
        ELSE
          D(1)=0.D0
          D(2)=0.D0
          D(3)=0.D0
          D(4)=0.D0
        ENDIF
      CLOSE(UNIT=20)


      CALL FIND_PICFILE(PICNM,PICTFILE,EX)
      IF(.NOT.EX) RETURN

      do i=-qsz,qsz
      do j=-qsz,qsz
        dnx(i,j)=-1
        if(huse(i,j)) then
          w(1)=v(1)+scale*(i*uy(1)+j*ux(1)+ht(i,j)*uz(1))
          w(2)=v(2)+scale*(i*uy(2)+j*ux(2)+ht(i,j)*uz(2))
          w(3)=v(3)+scale*(i*uy(3)+j*ux(3)+ht(i,j)*uz(3))
          CALL V2IMGPL(w,V0,PICNM,NPX,NLN,MMFL,CTR,KMAT,D,
     .                     CX,CY,CZ, EX,z)
          if(ex) dnx(i,j)=brt0(Z,npx,nln,dn)
        endif
      enddo
      enddo

C JRW comment. The below looks for the highest value.
      z2=0
      do i=-qsz,qsz
      do j=-qsz,qsz
      if(dnx(i,j).ge.0) then
        real4=z2
        z2=max(real4,dnx(i,j))
      endif
      enddo
      enddo

C JRW comment. The below scales to the max value. Commenting out.
C      do i=-qsz,qsz
C      do j=-qsz,qsz
C      if(dnx(i,j).gt.0) then
C        dnx(i,j)=dnx(i,j)/z2
C      endif
C      enddo
C      enddo

      return
      end

c  ...............................................................
      SUBROUTINE MM2VN(MM,MMFL,V0,CX,CY,CZ,USE,V,N)
c  ...............................................................

      IMPLICIT NONE

      real*8        w1, w2, a, b, c, d0,
     .              x(0:3), y(0:3), z(0:3), eps,
     .              v0dotux, v0dotuy, v0dotuz, p0, p1, p2, p3,  
     .              vec(3,0:512,0:512,6), a0(0:512,0:512,6)

      integer*4     q, i, j, f, l, i0, j0, f0, i1, j1, iff
      
      DOUBLE PRECISION      VDOT
      DOUBLE PRECISION      VNORM
      DOUBLE PRECISION      V0(3)
      DOUBLE PRECISION      V(3)
      DOUBLE PRECISION      N(3)
      DOUBLE PRECISION      CX(3)      
      DOUBLE PRECISION      CY(3)
      DOUBLE PRECISION      CZ(3)
      DOUBLE PRECISION      UX(3)      
      DOUBLE PRECISION      UY(3)
      DOUBLE PRECISION      UZ(3)
      DOUBLE PRECISION      MMFL
      DOUBLE PRECISION      MM(2)
      DOUBLE PRECISION      mx2

      CHARACTER*72          INFILE

      LOGICAL               USE

      save         iff, mx2, q, vec
      data         iff/0/

      if(iff.eq.0) then
        iff=1
        INFILE='SHAPEFILES/SHAPE.TXT'
        call get_model(infile,q,vec,a0)
        mx2=0
        do f=1,6
        do j=0,q-1
        do i=0,q-1
          CALL VSUB(VEC(1,I,J,F),VEC(1,I+1,J+1,F),V)
          mx2=MAX(mx2,VNORM(V))
          CALL VSUB(VEC(1,I+1,J,F),VEC(1,I,J+1,F),V)
          mx2=MAX(mx2,VNORM(V))
        enddo
        enddo
        enddo
        mx2=mx2**2
      endif

      DO I=1,3
        UZ(I)=MM(1)*CX(I)+MM(2)*CY(I)+MMFL*CZ(I)
      ENDDO
      CALL VHAT(UZ,UZ)
      CALL UCRSS(UZ,CX,UY)
      CALL UCRSS(UY,UZ,UX)

      v0dotux=VDOT(V0,UX)
      v0dotuy=VDOT(V0,UY)
      v0dotuz=VDOT(V0,UZ)

      d0=1.d10
      USE=.FALSE.
      eps=0.01

      do f=1,6
      do j=0,q-1
      do i=0,q-1
        l=-1
        do j1=0,1
        do i1=0,1
          l=l+1
          call vequ(vec(1,i+i1,j+j1,f),v)
          x(l)=VDOT(ux,v)+v0dotux
          y(l)=VDOT(uy,v)+v0dotuy
          z(l)=VDOT(uz,v)+v0dotuz
          IF((x(l)**2+y(l)**2).gt.mx2) GO TO 10
        enddo
        enddo
        p0=-x(0)*y(1)+y(0)*x(1)
        p1=-x(1)*y(3)+y(1)*x(3)
        p2=-x(3)*y(2)+y(3)*x(2)
        p3=-x(2)*y(0)+y(2)*x(0)
        if((p0.le.eps).and.(p1.le.eps).and.
     .     (p2.le.eps).and.(p3.le.eps)) then
          x(3)=x(0)-x(1)-x(2)+x(3)
          x(1)=x(1)-x(0)
          x(2)=x(2)-x(0)
          y(3)=y(0)-y(1)-y(2)+y(3)
          y(1)=y(1)-y(0)
          y(2)=y(2)-y(0)
          z(3)=z(0)-z(1)-z(2)+z(3)
          z(1)=z(1)-z(0)
          z(2)=z(2)-z(0)
          a=x(1)*y(3)-x(3)*y(1)
          b=x(1)*y(2)-x(2)*y(1)-y(0)*x(3)+x(0)*y(3)
          c=-y(0)*x(2)+x(0)*y(2)
          if(b**2.gt.4*a*c) then
            if(abs(a).gt.(1.d-8)) then
              w1=(sqrt(1-4*a*c/b**2)-1)*b/(2*a)
            else
              w1=-c/b-a*c**2/b**3
            endif
            if(abs(x(2)+x(3)*w1).gt.abs(y(2)+y(3)*w1)) then
              w2=(-x(0)-x(1)*w1)/(x(2)+x(3)*w1)
            else
              w2=(-y(0)-y(1)*w1)/(y(2)+y(3)*w1)
            endif
            if((z(0)+z(1)*w1+z(2)*w2+z(3)*w1*w2).lt.d0) then
              d0=z(0)+z(1)*w1+z(2)*w2+z(3)*w1*w2
              i0=i
              j0=j
              f0=f
              USE=.TRUE.
            endif
          endif
        endif
10      continue
      enddo
      enddo
      enddo

      i0=0
      j0=0
      f0=1
      IF(USE) THEN
        v(1)=D0*UZ(1)-V0(1)
        v(2)=D0*UZ(2)-V0(2)
        v(3)=D0*UZ(3)-V0(3)
        CALL VSUB(VEC(1,I0+1,J0+1,F0),VEC(1,I0,J0,F0),UX)
        CALL VSUB(VEC(1,I0+1,J0,F0),VEC(1,I0,J0+1,F0),UY)
        CALL UCRSS(UX,UY,N)
      ENDIF

      return
      end

C   ..................................................
      SUBROUTINE PXMM(PICNM,MMFL,KMAT,CTR,D,PX,MM)
C   ..................................................

      IMPLICIT NONE
          
      DOUBLE PRECISION      KMAT(2,3)
      DOUBLE PRECISION      CTR(2)
      DOUBLE PRECISION      MMFL
      DOUBLE PRECISION      Z(2)  
      DOUBLE PRECISION      MM(2)  
      DOUBLE PRECISION      PX(2)
      DOUBLE PRECISION      PX0(2)
      DOUBLE PRECISION      D(4)
      DOUBLE PRECISION      Z1, Z2, Z3
      INTEGER               I
      CHARACTER*12          PICNM
          
      MM(1)=0
      MM(2)=0
      PX0(1)=CTR(1)
      PX0(2)=CTR(2)
      
      DO I=1,20
        MM(1)=MM(1)+(PX(1)-PX0(1))/KMAT(1,1)
        MM(2)=MM(2)+(PX(2)-PX0(2))/KMAT(2,2)
        Z1=MM(1)
        Z2=MM(2)
        Z3=MMFL
        call getmm(picnm,mmfl,z1,z2,z3,z)
        CALL MMPX(KMAT,CTR,D,Z,PX0)
        IF(((PX(1)-PX0(1))**2
     .    + (PX(2)-PX0(2))**2).LT.(1.d-8)) RETURN
      ENDDO

      RETURN
      END

C   ..................................................
      FUNCTION SLEN(STRING)
C   ..................................................

      IMPLICIT NONE
      
      INTEGER*4        SLEN, I
      CHARACTER*(*)     STRING
      
      SLEN=LEN(STRING)
      DO I=1,SLEN
      IF(STRING(I:I).EQ.' ') THEN
        SLEN=I-1
        RETURN
      ENDIF
      ENDDO

      RETURN
      END            
      
c   ..................................................
      subroutine get_model(infile,q,vec,a0)
c   ..................................................

      IMPLICIT NONE

      REAL*8             VEC(3,0:512,0:512,6), A0(0:512,0:512,6)
      INTEGER*4          I, J, K, F, Q, N
      CHARACTER*72       INFILE
      CHARACTER*100      LINE

C     Check to see if file contain albedos, i.e. N=4
      N=0
      OPEN(UNIT=49,FILE=INFILE,STATUS='OLD')
      READ(49,*)
      READ(49,FMT='(A100)') LINE
      DO K=2,100
        IF((LINE(K-1:K-1).NE.' ').AND.(LINE(K:K).EQ.' ')) N=N+1
      ENDDO
      CLOSE(UNIT=49)

C     Re-open shape model file and read vector values
      OPEN(UNIT=49,FILE=INFILE,STATUS='OLD')
      READ(49,*) Q
      DO F=1,6
        DO J = 0,Q
          DO I = 0,Q
            A0(I,J,F)=1.D0
            READ(49,FMT='(A100)') LINE
            IF(N.EQ.3) THEN
              READ(LINE,*) (VEC(K,I,J,F), K=1,3)
            ENDIF
            IF(N.EQ.4) THEN
              READ(LINE,*) (VEC(K,I,J,F), K=1,3), A0(I,J,F)
            ENDIF
          ENDDO
        ENDDO
      ENDDO
      CLOSE(UNIT=49)
      
      RETURN
      END
      
c  ..............................................
      subroutine raw2pgm(infile,outfile,npx,nln)
c  ..............................................

      implicit none
      
      integer*4             npx, nln, i, j, k
      character*10          line
      character*50          header
      character*72          infile, outfile
      character*30000000    dn

      
      header(1:2)='P5'
      header(3:3)=char(10)
      header(4:5)='#.'
      header(6:6)=char(10)
      k=7
      write(line,fmt='(i10)') npx
      do i=1,10
      if(line(i:i).ne.' ') then
        header(k:k+10-i) = line(i:10)
        k=k+11-i
        go to 10
      endif
      enddo
10    continue
      header(k:k)=' '
      k=k+1
      write(line,fmt='(i10)') nln
      do i=1,10
      if(line(i:i).ne.' ') then
        header(k:k+10-i) = line(i:10)
        k=k+11-i
        go to 20
      endif
      enddo
20    continue
      header(k:k)=char(10)
      k=k+1
      header(k:k+2)='255'
      k=k+3
      header(k:k)=char(10)
      
      open(unit=10, file=infile, recl=npx, access='direct', 
     .     status='old')
        dn(1:k)=header(1:k)
        do j=1,nln
          read(10,rec=j) dn(k+(j-1)*npx+1:k+j*npx)
        enddo
      close(unit=10)
      open(unit=20, file=outfile, status='unknown')
      close(unit=20, status='delete')
      open(unit=20, file=outfile, recl=npx, access='direct', 
     .     status='unknown')
        do j=1,nln
          write(20,rec=j) dn(1+(j-1)*npx:j*npx)
        enddo
      close(unit=20)
      open(unit=20, file=outfile, recl=1, access='direct', 
     .     status='unknown')
        do j=npx*nln+1,npx*nln+k
          write(20,rec=j) dn(j:j)
        enddo
      close(unit=20)
      
      return
      end

c   ................................................
      SUBROUTINE READ_MAP(LMRKFILE,NTMP,QSZ,SCALE,V,UX,UY,UZ,HT,ALB)
c   ................................................

      IMPLICIT NONE
      
      INTEGER               NTMP
      INTEGER               QSZ
      INTEGER               I
      INTEGER               J
      INTEGER               K
      INTEGER               K0
      INTEGER               IX(24)
      INTEGER               JX(24)
      INTEGER               NREC

      DOUBLE PRECISION      SCALE
      DOUBLE PRECISION      HSCALE
      DOUBLE PRECISION      V(3)
      DOUBLE PRECISION      UX(3)
      DOUBLE PRECISION      UY(3)
      DOUBLE PRECISION      UZ(3)

      REAL*4                HT(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                ALB(-NTMP:NTMP,-NTMP:NTMP)

      CHARACTER*72          BLINE
      CHARACTER*72          LMRKFILE

      CHARACTER*1           CH1
      CHARACTER*2           CH2, CH2F
      CHARACTER*4           CH4, CH4F
      INTEGER*2             IX2
      REAL*4                RL4
      EQUIVALENCE          (IX2,CH2)
      EQUIVALENCE          (RL4,CH4)

      CHARACTER*2           C2
      INTEGER*2             I2
      EQUIVALENCE          (I2,C2)
      LOGICAL               LFLAG

      c2='69'
      LFLAG=.TRUE.
      if(i2.eq.13881) LFLAG=.FALSE.

      OPEN(UNIT=20, FILE=LMRKFILE, ACCESS='DIRECT',
     .     RECL=72, status='OLD')

        READ(20,REC=1) BLINE
        CH4f=BLINE(7:10)
        call flip(4,lflag,ch4f,ch4)
        SCALE=RL4
        QSZ=ICHAR(BLINE(11:11))
     .     +ICHAR(BLINE(12:12))*256
        DO K=1,3
          CH4f=BLINE(12+4*K:15+4*K)
          call flip(4,lflag,ch4f,ch4)
          V(K)=RL4
          CH4f=BLINE(24+4*K:27+4*K)
          call flip(4,lflag,ch4f,ch4)
          UX(K)=RL4
          CH4f=BLINE(36+4*K:39+4*K)
          call flip(4,lflag,ch4f,ch4)
          UY(K)=RL4
          CH4f=BLINE(48+4*K:51+4*K)
          call flip(4,lflag,ch4f,ch4)
          UZ(K)=RL4
        ENDDO
        CH4f=BLINE(64:67)
        call flip(4,lflag,ch4f,ch4)
        HSCALE=RL4

        DO I=-QSZ,QSZ
        DO J=-QSZ,QSZ
          ALB(I,J)=0
          HT(I,J)=0
        ENDDO
        ENDDO
        
        NREC=1        
        K=0
        DO J=-QSZ,QSZ
        DO I=-QSZ,QSZ
          K=K+1
          IX(K)=I
          JX(K)=J
          IF(K.EQ.24) THEN
            NREC=NREC+1
            READ(20,REC=NREC) BLINE
            DO K=1,24
              CH1=BLINE(3*K:3*K)
              IF(ICHAR(CH1).NE.0) THEN
                CH2f=BLINE(3*K-2:3*K-1)
                call flip(2,lflag,ch2f,ch2)
                HT(IX(K),JX(K))=HSCALE*IX2
                ALB(IX(K),JX(K))=.01*ICHAR(CH1)
              ENDIF
            ENDDO
            K=0
          ENDIF
        enddo
        enddo
        IF(K.NE.0) THEN
          K0=K
          NREC=NREC+1
          READ(20,REC=NREC) BLINE
          DO K=1,K0
            CH1=BLINE(3*K:3*K)
            CH2f=BLINE(3*K-2:3*K-1)
            call flip(2,lflag,ch2f,ch2)
            HT(IX(K),JX(K))=HSCALE*IX2
            ALB(IX(K),JX(K))=.01*ICHAR(CH1)
          ENDDO
        ENDIF

      CLOSE(UNIT=20)
      
      RETURN
      END     

c   ..................................................
      subroutine flip(n,lflag,ch1,ch2)
c   ..................................................

      integer*4        n, i
      character*(*)    ch1, ch2
      logical          lflag

      if(lflag) then
        do i=1,n
          ch2(i:i)=ch1(n-i+1:n-i+1)
        enddo
      else
        ch2=ch1
      endif
      
      return
      end
      

