c  g77 PROCESS_IMG.f /Users/rgaskell/toolkit/lib/spicelib.a -o PROCESS_IMG.e
c  ftn PROCESS_IMG.f /Users/rgaskell/toolkit/lib/spicelib.a -o PROCESS_IMG.e
c  gfortran processLROeep.f /usr/local/toolkit/lib/spicelib.a -o ~/bin/t.processLRO

C  Sent 10 October 2013.  Used to process LRO data.  
C		Also makes sumfiles.  The 32 bit images are made into 16 bit
C	9 Jan 2014 - renamed
C 17 July 2019 uncommented so DN < 0 becomes DN=0

      IMPLICIT NONE

      CHARACTER*(*)         D2COL
      PARAMETER           ( D2COL = '(2X,A11,E23.16,5X,E23.16)' )

      INTEGER               STRLEN
      PARAMETER           ( STRLEN =  32 )

      DOUBLE PRECISION      RPD
      DOUBLE PRECISION      SPD
      DOUBLE PRECISION      VDOT
      DOUBLE PRECISION      VNORM
      DOUBLE PRECISION      EXP
      DOUBLE PRECISION      IFD
      DOUBLE PRECISION      LT
      DOUBLE PRECISION      ET, ET0, ET1, ETL
      DOUBLE PRECISION      ETX0, ETX1
      DOUBLE PRECISION      SCLK0, SCLK1
      DOUBLE PRECISION      X1, Y1, Z1
      DOUBLE PRECISION      X2, Y2, Z2
      DOUBLE PRECISION      X3, Y3, Z3
      DOUBLE PRECISION      a(-1:1)
      DOUBLE PRECISION      b
      DOUBLE PRECISION      c
      DOUBLE PRECISION      dpx_dln(60000)
      DOUBLE PRECISION      dpx(60000)
      DOUBLE PRECISION      KMAT(2,3)
      DOUBLE PRECISION      CPIC(2)
      DOUBLE PRECISION      MMFL
      DOUBLE PRECISION      K11
      DOUBLE PRECISION      D(4)
      DOUBLE PRECISION      SIG_VSO(3)
      DOUBLE PRECISION      SIG_PTG(3)
      DOUBLE PRECISION      X(2000)
      DOUBLE PRECISION      Y(2000)
      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      SEC_PER_LINE
      DOUBLE PRECISION      MAT(20,20)
      DOUBLE PRECISION      IMAT(20,20)
      DOUBLE PRECISION      VDUM(6)
      DOUBLE PRECISION      USUN(3)
      DOUBLE PRECISION      U(3)
      DOUBLE PRECISION      W(3)
      DOUBLE PRECISION      V(3)
      DOUBLE PRECISION      V0(3)
      DOUBLE PRECISION      VX(3)
      DOUBLE PRECISION      VY(3)
      DOUBLE PRECISION      VZ(3)
      DOUBLE PRECISION      WX(3)
      DOUBLE PRECISION      WY(3)
      DOUBLE PRECISION      WZ(3)
      DOUBLE PRECISION      CX(3)
      DOUBLE PRECISION      CY(3)
      DOUBLE PRECISION      CZ(3)
      DOUBLE PRECISION      CPX
      DOUBLE PRECISION      WMAT(20,3)    
      DOUBLE PRECISION      CMAT(20,3)    
      DOUBLE PRECISION      CAX(20,3)
      DOUBLE PRECISION      MA(3,3)
      DOUBLE PRECISION      MB(3,3)
      DOUBLE PRECISION      R00
      DOUBLE PRECISION      WW(8)
      DOUBLE PRECISION      XX(8)
      DOUBLE PRECISION      PAR(4)
      DOUBLE PRECISION      DSCOBJ
      DOUBLE PRECISION      DPTG(3)
      DOUBLE PRECISION      ZW(7)
      DOUBLE PRECISION      ROTP


      INTEGER               OBSRVR
      INTEGER               TARGET
      INTEGER               SLEN
      INTEGER               NLF
      INTEGER               NEQ
      INTEGER               I
      INTEGER               IMAX
      INTEGER               J
      INTEGER               K
      INTEGER               K0
      INTEGER               K1
      INTEGER               K2
      INTEGER               KMAX
      INTEGER               L
      INTEGER               N
      INTEGER               NC
      INTEGER               NS
      INTEGER               SCL
      INTEGER               NPX
      INTEGER               NLN
      INTEGER               NLN0
      INTEGER               NPX1
      INTEGER               NLN1
      INTEGER               RECL
      INTEGER               BITS
      INTEGER               DN(6000,60000)
      INTEGER               DN1(2000,20000)
      INTEGER               DN2(2000,20000)
      INTEGER               PREC
      INTEGER               RECORD_BYTES
      INTEGER               FILE_RECORDS
      INTEGER               LABEL_RECORDS
      INTEGER               NFLT
      INTEGER               OFFSET
      INTEGER               NFM
      INTEGER               KF
      INTEGER               HGT(0:7)
      INTEGER               CYCLE
      INTEGER               T1, T2
      INTEGER               DNMAX, DNMIN

      CHARACTER*1           ANS
      CHARACTER*1           SUF
      CHARACTER*1           NOSCLK
      CHARACTER*3           FORMAT
      CHARACTER*3           FLT(7)
      CHARACTER*4           VERSION
      CHARACTER*10          TM
      CHARACTER*12          PICNM
      CHARACTER*13          XNAME(100000)
      CHARACTER*17          IMGNM
      CHARACTER*24          UTC
      CHARACTER*24          UTC0
      CHARACTER*24          UTC1
      CHARACTER*17          CLK
      CHARACTER*72          PICFILE
      CHARACTER*72          INFILE
      CHARACTER*72          TMPFILE
      CHARACTER*80          ZLINE0
      CHARACTER*80          ZLINE(100000)
      CHARACTER*80          KERNEL
      CHARACTER*50000       LINE
      CHARACTER*1024        YLINE
      CHARACTER*2048        XLINE
      CHARACTER*(STRLEN)    INST
      CHARACTER*(STRLEN)    BFRAME
      CHARACTER*(STRLEN)    ABCORR

      LOGICAL               EX
      LOGICAL               FOUND(2000)
      LOGICAL               ONPX
      LOGICAL               FLIP
      LOGICAL               EXTRA_FLIP

      DOUBLE PRECISION      DEG_PER_DAY
      PARAMETER            (DEG_PER_DAY=13.1763581)

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

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

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

      CALL RESLUN(10)
      CALL RESLUN(20)
      CALL RESLUN(22)
      CALL RESLUN(25)
      CALL RESLUN(69)

      FORMAT = 'C'
      PREC   =  3
      ABCORR='LT+S'
      TARGET=301 
      BFRAME='MOON_ME'
      OBSRVR=-85
      POL_AXIS(1)=0.D0
      POL_AXIS(2)=0.D0
      POL_AXIS(3)=1.D0
      T1=5
      T2=255
      DO I=1,2
      DO J=1,3
        KMAT(I,J)=0
      ENDDO
      ENDDO
      DO I=1,4
        D(I)=0
      ENDDO  
      DO I=1,3
        SIG_VSO(I)=1.D-3
        SIG_PTG(I)=2.D-4
      ENDDO
      SIG_PTG(3)=1.D-5
      R00=1737.4

      WRITE(6,*) 'LOADING KERNELS'
      OPEN(UNIT=10, FILE='../LRO/PROCESS_IMG.IN', STATUS='OLD')
10      CONTINUE
        READ(10,FMT='(A80)') KERNEL
        IF(KERNEL(1:1).EQ.'#') GO TO 10
        IF(KERNEL(1:2).EQ.'ZW') THEN
          READ(KERNEL(3:80),*) K, Z1
          ZW(K)=Z1
          GO TO 10
        ENDIF
        INQUIRE(FILE=KERNEL, EXIST=EX)
        IF(EX) THEN
          CALL FURNSH ( KERNEL )
        ELSE 
          write (*,*) 'Cant find: ', KERNEL
        ENDIF
        IF(KERNEL(1:11).NE.'END_KERNELS') GO TO 10
      CLOSE(UNIT=10)

30    CONTINUE
      WRITE(6,*) 
      WRITE(6,*) 'Input imgnm'
      READ(5,FMT='(A17)') IMGNM
      PICNM=IMGNM(1:11)//' '

      INFILE='NEW_IMAGES/'//IMGNM
      write (*,*) infile
      write (*,*) picnm
      INQUIRE(FILE=INFILE, EXIST=EX)
      IF(.NOT.EX) THEN
        WRITE(6,*) 'no image file'
        STOP
      ENDIF

      NOSCLK='n'
C      WRITE(6,*) 'Veto sclk? (y/n)'
C      READ(5,FMT='(A1)') NOSCLK

      SCL=1
      NFM=0
      OPEN(UNIT=10,FILE=INFILE,ACCESS='DIRECT',RECL=1000, STATUS='OLD')
        READ(10,REC=1) LINE(1:1000)
        DO I=1,1000
          IF(LINE(I:I+11).EQ.'RECORD_BYTES') THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            READ(LINE(J:K),*) RECORD_BYTES
          ENDIF
          IF(LINE(I:I+11).EQ.'FILE_RECORDS') THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            READ(LINE(J:K),*) FILE_RECORDS
          ENDIF
          IF(LINE(I:I+12).EQ.'LABEL_RECORDS') THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            READ(LINE(J:K),*) LABEL_RECORDS
          ENDIF
        ENDDO
      CLOSE(UNIT=10)

      RECL=LABEL_RECORDS*RECORD_BYTES
      OPEN(UNIT=10,FILE=INFILE,ACCESS='DIRECT',RECL=RECL,STATUS='OLD')
      write (*,*) "Open ", infile
        READ(10,REC=1) LINE(1:RECL)
        DO I=1,RECL
          IF(NOSCLK.NE.'y') THEN
            IF(LINE(I:I+27).EQ.'SPACECRAFT_CLOCK_START_COUNT') THEN
              J=NEQ(LINE,I)+1
              K=NLF(LINE,I)-1
              IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
              K0=0
              DO K1=J,K
                IF(LINE(K1:K1).EQ.'/') K0=K1
              ENDDO
              IF(K0.NE.0) THEN
                CLK=LINE(K0-1:K0+15)
                CALL SCENCD (OBSRVR, CLK, SCLK0 )
              ELSE
                READ(LINE(J:K),*) SCLK0
              ENDIF
              CALL SCT2E (OBSRVR, SCLK0, ET0 )
            ENDIF
            IF(LINE(I:I+26).EQ.'SPACECRAFT_CLOCK_STOP_COUNT') THEN
              J=NEQ(LINE,I)+1
              K=NLF(LINE,I)-1
              IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
              K0=0
              DO K1=J,K
                IF(LINE(K1:K1).EQ.'/') K0=K1
              ENDDO
              IF(K0.NE.0) THEN
                CLK=LINE(K0-1:K0+15)
                CALL SCENCD (OBSRVR, CLK, SCLK1 )
              ELSE
                READ(LINE(J:K),*) SCLK1
              ENDIF
              CALL SCT2E (OBSRVR, SCLK1, ET1 )
            ENDIF
          ENDIF
          IF(LINE(I:I+9).EQ.'START_TIME') THEN
            J=NEQ(LINE,I)+1
            UTC0=LINE(J:J+23)
            CALL UTC2ET ( UTC0, ETX0 )
          ENDIF
          IF(LINE(I:I+8).EQ.'STOP_TIME') THEN
            J=NEQ(LINE,I)+1
            UTC1=LINE(J:J+23)
            CALL UTC2ET ( UTC1, ETX1 )
          ENDIF
          IF(LINE(I:I+17).EQ.'CROSSTRACK_SUMMING') THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            READ(LINE(J:K),*) SCL
          ENDIF
          IF(LINE(I:I+4).EQ.'LINES') THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            READ(LINE(J:K),*) NLN0
            NLN=NLN0
          ENDIF
          IF(LINE(I:I+11).EQ.'LINE_SAMPLES') THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            READ(LINE(J:K),*) NPX
          ENDIF
          IF(LINE(I:I+10).EQ.'SAMPLE_BITS') THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            READ(LINE(J:K),*) BITS
          ENDIF
          IF(LINE(I:I+10).EQ.'LRO:NFRAMES') THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            READ(LINE(J:K),*) NFM
          ENDIF
          IF(LINE(I:I+16).EQ.'EXPOSURE_DURATION') THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            READ(LINE(J:K),*) EXP
          ENDIF
          IF(LINE(I:I+15).EQ.'INTERFRAME_DELAY') THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            READ(LINE(J:K),*) IFD
          ENDIF
          IF(LINE(I:I+10).EQ.'FILTER_NAME') THEN
          IF((PICNM(11:11).EQ.'C').OR.(PICNM(11:11).EQ.'M')) THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            K0=0
            DO K1=J,K
            IF((ICHAR(LINE(K1:K1)).GE.48).AND.
     .         (ICHAR(LINE(K1:K1)).LE.57).AND.
     .         (ICHAR(LINE(K1+1:K1+1)).GE.48).AND.
     .         (ICHAR(LINE(K1+1:K1+1)).LE.57).AND.
     .         (ICHAR(LINE(K1+2:K1+2)).GE.48).AND.
     .         (ICHAR(LINE(K1+2:K1+2)).LE.57)) THEN
              K0=K0+1
              FLT(K0)=LINE(K1:K1+2)
              HGT(K0)=14
              IF(FLT(K0).EQ.'320') FLT(K0)='321'
              IF(FLT(K0).EQ.'565') FLT(K0)='566'
              IF(FLT(K0).EQ.'605') FLT(K0)='604'
              IF(FLT(K0).EQ.'645') FLT(K0)='643'
              IF(FLT(K0).EQ.'690') FLT(K0)='689'
              IF((FLT(K0).EQ.'321').OR.(FLT(K0).EQ.'360')) HGT(K0)=4
            ENDIF
            ENDDO
            NFLT=K0
          ENDIF
          ENDIF
          IF(LINE(I:I+23).EQ.'CENTER_FILTER_WAVELENGTH') THEN
          IF((PICNM(11:11).EQ.'C').OR.(PICNM(11:11).EQ.'M')) THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            K0=0
            DO K1=J,K
            IF((ICHAR(LINE(K1:K1)).GE.48).AND.
     .         (ICHAR(LINE(K1:K1)).LE.57).AND.
     .         (ICHAR(LINE(K1+1:K1+1)).GE.48).AND.
     .         (ICHAR(LINE(K1+1:K1+1)).LE.57).AND.
     .         (ICHAR(LINE(K1+2:K1+2)).GE.48).AND.
     .         (ICHAR(LINE(K1+2:K1+2)).LE.57)) THEN
              K0=K0+1
              FLT(K0)=LINE(K1:K1+2)
              HGT(K0)=14
              IF((FLT(K0).EQ.'321').OR.(FLT(K0).EQ.'360')) HGT(K0)=4
            ENDIF
            ENDDO
            NFLT=K0
          ENDIF
          ENDIF
          IF(LINE(I:I+18).EQ.'PRODUCT_VERSION_ID') THEN
            J=NEQ(LINE,I)+1
            K=NLF(LINE,I)-1
            IF(ICHAR(LINE(K:K)).EQ.13) K=K-1
            IF(J.EQ.K) THEN
              VERSION='v0.0'
            ELSE
              READ(LINE(J:K),*) VERSION
            ENDIF
          ENDIF
        ENDDO
      CLOSE(UNIT=10)

      EXTRA_FLIP=(VERSION.GE.'v1.3').OR.(BITS.GT.8)

      IF(PICNM(11:11).EQ.'R') PICNM(1:1)='N'
      IF(PICNM(11:11).EQ.'L') PICNM(1:1)='N'
      IF(PICNM(11:11).EQ.'C') PICNM(1:1)='W'
      IF(PICNM(11:11).EQ.'M') PICNM(1:1)='W'

      IF(PICNM(11:11).EQ.'R') INST='LRO_LROCNACR'
      IF(PICNM(11:11).EQ.'L') INST='LRO_LROCNACL'
      IF(PICNM(1:1).EQ.'W')   INST='LRO_LROCWAC'

      IF((PICNM(1:1).EQ.'W').AND.(NFM.EQ.0)) THEN
        CYCLE=0
        DO K=1,NFLT
          CYCLE=CYCLE+HGT(K)
        ENDDO
        NFM=NLN/CYCLE
      ENDIF

      IF(PICNM(1:1).EQ.'W') THEN
      write (*,*) "Status W "
        K=(NFM+1)/2
        ET=((NFM-K)*ETX0+(K-1)*ETX1)/(NFM-1)
        SIG_PTG(1)=2.D-2                            !!!
      ELSE
        ET=(ETX0+ETX1)/2
      ENDIF
      CALL ET2UTC ( ET, FORMAT, PREC, UTC )
      CALL ET2UTC ( ETX0, FORMAT, PREC, UTC0 )
      CALL ET2UTC ( ETX1, FORMAT, PREC, UTC1 )
      write(6,*) 
      write(6,*) 'Times from onboard utc'
      write(6,*) 'UTC0 = ', utc0
      write(6,*) 'UTC  = ', utc
      write(6,*) 'UTC1 = ', utc1
      write(6,*) 'Frame time = ', ETX1-ETX0
      IF(NOSCLK.EQ.'y') THEN
        ET0=ETX0
        ET1=ETX1
      ELSE
        IF(PICNM(1:1).EQ.'W') THEN
          K=(NFM+1)/2
          ET=((NFM-K)*ET0+(K-1)*ET1)/(NFM-1)
        ELSE
          ET=(ET0+ET1)/2
        ENDIF
        CALL ET2UTC ( ET, FORMAT, PREC, UTC )
        CALL ET2UTC ( ET0, FORMAT, PREC, UTC0 )
        CALL ET2UTC ( ET1, FORMAT, PREC, UTC1 )
        write(6,*) 
        write(6,*) 'Times from spacecraft clock'
        write(6,*) 'UTC0 = ', utc0
        write(6,*) 'UTC  = ', utc
        write(6,*) 'UTC1 = ', utc1
        write(6,*) 'Frame time = ', ET1-ET0
        write(6,*) 
      ENDIF

      CALL SPKEZ (TARGET, ET, BFRAME, ABCORR, OBSRVR, VDUM, LT)
      CALL VEQU(VDUM(4),V)
      CALL VHAT(V,V)
      CALL PXFORM (INST, BFRAME, ET, MA)
      DO I=1,3
        CY(I)=-MA(I,1)
      ENDDO
      IF(VDOT(CY,V).GT.0) THEN
        FLIP=.TRUE.
        WRITE(6,*) 'FLIP'
      ELSE
        FLIP=.FALSE.
        WRITE(6,*) 'NO FLIP'
      ENDIF
      WRITE(6,*)

c      go to 30  ! Uncomment to check images for flip only.

      IF(PICNM(11:11).EQ.'C') THEN
      write (*,*) "Status C "
        DO K=1,NFLT
          WRITE(6,FMT='(I5, 5X, A3)') K, FLT(K)
        ENDDO
        WRITE(6,*) 'Choose filter'
        READ(5,*) KF
      write (*,*) "filter ", KF
        IF(FLT(KF).EQ.'321') SUF='1'
        IF(FLT(KF).EQ.'360') SUF='2'
        IF(FLT(KF).EQ.'415') SUF='3'
        IF(FLT(KF).EQ.'566') SUF='4'
        IF(FLT(KF).EQ.'604') SUF='5'
        IF(FLT(KF).EQ.'643') SUF='6'
        IF(FLT(KF).EQ.'689') SUF='7'
        OFFSET=0
        DO K=1,KF-1
          OFFSET=OFFSET+HGT(K)
        ENDDO
        CYCLE=0
        DO K=1,NFLT
          CYCLE=CYCLE+HGT(K)
        ENDDO
      ENDIF
      IF(PICNM(11:11).EQ.'M') THEN
        KF=1
        CYCLE=HGT(KF)
        OFFSET=0
        IF(FLT(KF).EQ.'321') SUF='1'
        IF(FLT(KF).EQ.'360') SUF='2'
        IF(FLT(KF).EQ.'415') SUF='3'
        IF(FLT(KF).EQ.'566') SUF='4'
        IF(FLT(KF).EQ.'604') SUF='5'
        IF(FLT(KF).EQ.'643') SUF='6'
        IF(FLT(KF).EQ.'689') SUF='7'
      ENDIF

35    CONTINUE
      RECL=RECORD_BYTES
      OPEN(UNIT=10,FILE=INFILE,ACCESS='DIRECT',RECL=RECL,
     .     STATUS='OLD')
      write (*,*) "Open ", infile, bits, flip, picnm(11:11)
        DO L=1,NLN
          IF(BITS.EQ.8) THEN
            READ(10,REC=L+LABEL_RECORDS) LINE(1:RECL)
          ENDIF
          IF(BITS.EQ.16) THEN
            READ(10,REC=2*L-1+LABEL_RECORDS) LINE(1:RECL)
            READ(10,REC=2*L+LABEL_RECORDS) LINE(RECL+1:2*RECL)
          ENDIF
          IF(BITS.EQ.32) THEN
            READ(10,REC=4*L-3+LABEL_RECORDS) LINE(1:RECL)
            READ(10,REC=4*L-2+LABEL_RECORDS) LINE(RECL+1:2*RECL)
            READ(10,REC=4*L-1+LABEL_RECORDS) LINE(2*RECL+1:3*RECL)
            READ(10,REC=4*L+LABEL_RECORDS) LINE(3*RECL+1:4*RECL)
          ENDIF
          DO I=1,NPX
            IF(BITS.EQ.8) THEN
              IX1=ICHAR(LINE(I:I))  
            ENDIF
            IF(BITS.EQ.16) THEN
              CH2F=LINE(2*I-1:2*I)
              CALL FLP(2,LSB,CH2F,CH2)
            ENDIF
            IF(BITS.EQ.32) THEN
              CH4F=LINE(4*I-3:4*I)
              CALL FLP(4,LSB,CH4F,CH4)
            ENDIF
            IF(FLIP) THEN
              IF(PICNM(11:11).EQ.'L') THEN 
                IF(BITS.EQ.8) DN(NPX+1-I,L)=IX1
                IF(BITS.EQ.16) DN(NPX+1-I,L)=IX2
              ENDIF
              IF((PICNM(11:11).EQ.'R').AND.(.NOT.EXTRA_FLIP)) THEN 
                IF(BITS.EQ.8) DN(NPX+1-I,L)=IX1
                IF(BITS.EQ.16) DN(NPX+1-I,L)=IX2
              ENDIF
              IF((PICNM(11:11).EQ.'R').AND.(EXTRA_FLIP)) THEN 
                IF(BITS.EQ.8) DN(I,L)=IX1
                IF(BITS.EQ.16) DN(I,L)=IX2
              ENDIF
              IF((PICNM(11:11).EQ.'C')) THEN
                K0=(L-1)/CYCLE
                K1=L-K0*CYCLE-OFFSET
                J=(K0+1)*HGT(KF)-K1+1
                IF((K1.GE.1).AND.(K1.LE.HGT(KF))) THEN
                  IF(BITS.EQ.8) DN(NPX+1-I,J)=IX1
                  IF(BITS.EQ.16) DN(NPX+1-I,J)=IX2
                  IF(BITS.EQ.32) DNX(NPX+1-I,J)=RL4
                ENDIF
              ENDIF
              IF((PICNM(11:11).EQ.'M')) THEN
                K0=(L-1)/CYCLE
                K1=L-K0*CYCLE
                J=(K0+1)*CYCLE-K1+1
                IF(BITS.EQ.8) DN(NPX+1-I,J)=IX1
                IF(BITS.EQ.16) DN(NPX+1-I,J)=IX2
                IF(BITS.EQ.32) DNX(NPX+1-I,J)=RL4
              ENDIF
            ELSE
              IF(PICNM(11:11).EQ.'L') THEN 
                IF(BITS.EQ.8) DN(I,L)=IX1
                IF(BITS.EQ.16) DN(I,L)=IX2
              ENDIF
              IF((PICNM(11:11).EQ.'R').AND.(.NOT.EXTRA_FLIP)) THEN 
                IF(BITS.EQ.8) DN(I,L)=IX1
                IF(BITS.EQ.16) DN(I,L)=IX2
              ENDIF
              IF((PICNM(11:11).EQ.'R').AND.(EXTRA_FLIP)) THEN 
                IF(BITS.EQ.8) DN(NPX+1-I,L)=IX1
                IF(BITS.EQ.16) DN(NPX+1-I,L)=IX2
              ENDIF
              IF((PICNM(11:11).EQ.'C')) THEN
                K0=(L-1)/CYCLE
                K1=L-K0*CYCLE-OFFSET
                J=K0*HGT(KF)+K1
                IF((K1.GE.1).AND.(K1.LE.HGT(KF))) THEN
                  IF(BITS.EQ.8)  DN(I,J)=IX1
                  IF(BITS.EQ.16) DN(I,J)=IX2
                  IF(BITS.EQ.32) DNX(I,J)=RL4
                ENDIF
              ENDIF
              IF((PICNM(11:11).EQ.'M')) THEN
                K0=(L-1)/CYCLE
                K1=L-K0*CYCLE
                J=K0*CYCLE+K1
                IF(BITS.EQ.8)  DN(I,J)=IX1
                IF(BITS.EQ.16) DN(I,J)=IX2
                IF(BITS.EQ.32) DNX(I,J)=RL4
              ENDIF
            ENDIF
          ENDDO
        ENDDO
      CLOSE(UNIT=10)
      IF((PICNM(11:11).EQ.'C')) THEN
        K0=(NLN-1)/CYCLE
        NLN=(K0+1)*HGT(KF)
      ENDIF

      IF(BITS.EQ.16) THEN
      write (*,*) "Bits == 16, commented"
        Z2=0
        DO J=1,NLN
        DO I=1,NPX
          IF(DN(I,J).LE.0) DN(I,J)=0
          Z2=MAX(DFLOAT(DN(I,J)),Z2)
        ENDDO
        ENDDO
        DO J=1,NLN
        DO I=1,NPX
c          DN(I,J)=NINT(32767*DN(I,J)/Z2)
        ENDDO
        ENDDO
        T1=640
        T2=32767
      ENDIF
      IF(BITS.EQ.32) THEN
        Z2=0
        DO J=1,NLN
        DO I=1,NPX
          IF(DNX(I,J).LT.0) DNX(I,J)=0
          Z2=MAX(DNX(I,J),Z2)
        ENDDO
        ENDDO
        DO J=1,NLN
        DO I=1,NPX
          DN(I,J)=NINT(32767*DNX(I,J)/Z2)
        ENDDO
        ENDDO
        T1=640
        T2=32767
      ENDIF

      IF(PICNM(11:11).EQ.'R') THEN
      write (*,*) "Status R "
        DO I=1,42/SCL
        DO J=1,NLN
          DN(NPX+1-I,J)=0
        ENDDO
        ENDDO
        DO I=NPX+1-24/SCL,NPX
        DO J=1,NLN
          DN(NPX+1-I,J)=0
        ENDDO
        ENDDO
        CPIC(1)=(2568+1)/SCL
        IF(FLIP) CPIC(1)=NPX+1-CPIC(1)
        MMFL=701.57        !NAC2     
        K11=142.8571/SCL
      ENDIF

      IF(PICNM(11:11).EQ.'L') THEN
      write (*,*) "Status L "
        DO I=1,42/SCL
        DO J=1,NLN
          DN(I,J)=0
        ENDDO
        ENDDO
        DO I=NPX+1-24/SCL,NPX
        DO J=1,NLN
          DN(I,J)=0
        ENDDO
        ENDDO
        CPIC(1)=(2548+1)/SCL
        IF(FLIP) CPIC(1)=NPX+1-CPIC(1)
        MMFL=699.62        !NAC1
        K11=142.8570/SCL
      ENDIF

      IF((PICNM(11:11).EQ.'C').AND.(SUF.LE.'2')) THEN
      write (*,*) "Status C 2 "
        NPX=128
        SCL=4
        DO I=1,NPX
        DO J=1,NLN
          DN(I,J)=DN(I+288,J)
        ENDDO
        ENDDO
        IF(FLIP) THEN
          DO J=1,NLN
            DN(1,J)=0
            DN(2,J)=0
          ENDDO
        ELSE
          DO J=1,NLN
            DN(NPX,J)=0
            DN(NPX-1,J)=0
          ENDDO
        ENDIF      
        PICNM(12:12)=SUF
c        MMFL=4.668         !WAC1
c        K11=111.055/SCL    !WAC1
        MMFL=4.693         !WAC2
        K11=111.097/SCL    !WAC2 
        CPIC(1)=254.3/SCL
        IF(FLIP) CPIC(1)=NPX+1-CPIC(1)
      ENDIF

      IF((PICNM(11:11).EQ.'C').AND.(SUF.GT.'2')) THEN
      write (*,*) "Status C  > 2"
        IF(FLIP) THEN
          DO J=1,NLN
            DN(NPX,J)=0
            DN(NPX-1,J)=0
          ENDDO
        ELSE
          DO J=1,NLN
            DN(1,J)=0
            DN(2,J)=0
          ENDDO
        ENDIF      
        PICNM(12:12)=SUF
c        MMFL=6.008         !WAC1
c        K11=111.111/SCL    !WAC1
        MMFL=6.013         !WAC2
        K11=111.019/SCL    !WAC2    
        CPIC(1)=349.9/SCL
        IF(FLIP) CPIC(1)=NPX+1-CPIC(1)
      ENDIF

      IF(PICNM(11:11).EQ.'M') THEN
      write (*,*) "Status M "
        IF(FLIP) THEN
          DO I=1,6/SCL
          DO J=1,NLN
            DN(I,J)=0
          ENDDO
          ENDDO
          DO I=NPX+1-10/SCL,NPX
          DO J=1,NLN
            DN(I,J)=0
          ENDDO
          ENDDO
        ELSE
          DO I=1,10/SCL
          DO J=1,NLN
            DN(I,J)=0
          ENDDO
          ENDDO
          DO I=NPX+1-6/SCL,NPX
          DO J=1,NLN
            DN(I,J)=0
          ENDDO
          ENDDO
        ENDIF
        PICNM(12:12)=SUF
c        MMFL=6.008         !WAC1
c        K11=111.111/SCL    !WAC1
        MMFL=6.013         !WAC2
        K11=111.019/SCL    !WAC2    
        CPIC(1)=509.9/SCL
        IF(FLIP) CPIC(1)=NPX+1-CPIC(1)
      ENDIF

      I=SLEN(PICNM)
      PICFILE='IMAGEFILES/'//PICNM(1:I)//'.DAT'
      IF(BITS.EQ.8) THEN
        OPEN(UNIT=10,FILE=PICFILE,ACCESS='DIRECT',RECL=NPX,
     .       STATUS='UNKNOWN')
      write (*,*) "write bits 8"
          DO J=1,NLN
            DO I=1,NPX
              LINE(I:I)=CHAR(DN(I,J))
            ENDDO
            WRITE(10,REC=J) LINE(1:NPX)
          ENDDO
        CLOSE(UNIT=10)
      ELSE
        OPEN(UNIT=10,FILE=PICFILE,ACCESS='DIRECT',RECL=2*NPX,
     .       STATUS='UNKNOWN')
      write (*,*) "write bits 16"
          DO J=1,NLN
            DO I=1,NPX
              K1=DN(I,J)/256
              K0=DN(I,J)-256*K1
              LINE(2*I-1:2*I-1)=CHAR(K0)
              LINE(2*I:2*I)=CHAR(K1)
            ENDDO
            WRITE(10,REC=J) LINE(1:2*NPX)
          ENDDO
        CLOSE(UNIT=10)
      ENDIF

      DNMIN=1000000
      DNMAX=0
      DO J=1,NLN
      DO I=1,NPX
        IF(DN(I,J).NE.0) THEN
          DNMAX=MAX(DNMAX,DN(I,J))
          DNMIN=MIN(DNMIN,DN(I,J))
        ENDIF
      ENDDO
      ENDDO

      OPEN(UNIT=10, FILE='TEMPFILE.gray',RECL=NPX,
     .     ACCESS='DIRECT', STATUS='UNKNOWN')
      write (*,*) "gray"
        DO J=1,NLN
          DO I=1,NPX
            LINE(I:I)=CHAR(0)
            IF(DN(I,J).NE.0) THEN
              LINE(I:I)=CHAR(NINT(1+254.*(DN(I,J)-DNMIN)/(DNMAX-DNMIN)))
            ENDIF
          ENDDO
          WRITE(10,REC=J) LINE(1:NPX)
        ENDDO
      CLOSE(UNIT=10)
      infile= 'TEMPFILE.gray' 
      tmpfile= 'TEMPFILE0.pgm'  
      call raw2pgm(infile,tmpfile,NPX,NLN)
      open(unit=63, file=infile, status='old')
      close(unit=63, status='delete')
      write(6,*) 'gc TEMPFILE0.pgm'
 
C      IF((PICNM(11:11).EQ.'C').OR.(PICNM(11:11).EQ.'M')) THEN
C        WRITE(6,*) 'Flip? (y/n)'
C        READ(5,FMT='(A1)') ANS
C        IF(ANS.EQ.'y') THEN
C          FLIP=.NOT.FLIP
C          NLN=NLN0
C          GO TO 35
C        ENDIF
C      ENDIF

      IF((PICNM(11:11).EQ.'L').OR.(PICNM(11:11).EQ.'R')) THEN
      write (*,*) "Status L  or R"
        CPIC(2)=(1+NLN)/2.D0
        SEC_PER_LINE=(ET1-ET0)/NLN
        WRITE(6,*) 'Jitter Fitter? (y/n)'
        READ(5,FMT='(A1)') ANS
        IF(ANS.NE.'y') NPX1=NPX
        IF(ANS.EQ.'y') THEN
          I=SLEN(PICNM)
          PICFILE='IMAGEFILES/'//PICNM(1:I)//'.DAT'
          OPEN(UNIT=10,FILE=PICFILE,STATUS='OLD')
          CLOSE(UNIT=10,STATUS='DELETE')
          PICNM(12:12)='X'
          do j=2,nln-1
            do k=-1,1
              a(k)=0.d0
            enddo
            do i=50,npx-50
            do k=-1,1
              z1=dn(i+k,j+1)+dn(i-k,j-1)-2*dn(i,j)
              a(k)=a(k)+z1**2
            enddo
            enddo
            b=(a(1)-a(-1))/2
            c=a(1)+a(-1)-2*a(0)
            dpx_dln(j)=b/(2*c)
          enddo
          dpx(1)=0.d0
          dpx_dln(1)=2*dpx_dln(2)-dpx_dln(3)        
          dpx_dln(nln)=2*dpx_dln(nln-1)-dpx_dln(nln-2)        
          z1=dpx(1)
          z2=dpx(1)
          I=SLEN(PICNM)
          do j=2,nln
            dpx(j)=dpx(j-1)+dpx_dln(j)
            z1=min(z1,dpx(j))
            z2=max(z2,dpx(j))
          enddo
          NPX1=NPX-NINT(Z1)+NINT(Z2)
          I=SLEN(PICNM)
          PICFILE='IMAGEFILES/'//PICNM(1:I)//'.DAT'
          IF(BITS.EQ.8) THEN
            OPEN(UNIT=10,FILE=PICFILE,ACCESS='DIRECT',RECL=NPX1,
     .           STATUS='UNKNOWN')
              DO J=1,NLN
                DO I=1,NPX1
                  LINE(I:I)=CHAR(0)
                ENDDO
                DO I=1,NPX
                  K=I-NINT(Z1)+NINT(DPX(J))
                  LINE(K:K)=CHAR(DN(I,J))
                ENDDO
                WRITE(10,REC=J) LINE(1:NPX1)
                DO I=1,NPX1
                  DN(I,J)=ICHAR(LINE(I:I))
                ENDDO
              ENDDO
            CLOSE(UNIT=10)
          ELSE
            OPEN(UNIT=10,FILE=PICFILE,ACCESS='DIRECT',RECL=2*NPX1,
     .           STATUS='UNKNOWN')
              DO J=1,NLN
                DO I=1,2*NPX1
                  LINE(I:I)=CHAR(0)
                ENDDO
                DO I=1,NPX
                  K=I-NINT(Z1)+NINT(DPX(J)) 
                  K1=DN(I,J)/256
                  K0=DN(I,J)-256*K1
                  LINE(2*K-1:2*K-1)=CHAR(K0)
                  LINE(2*K:2*K)=CHAR(K1)
                ENDDO
                WRITE(10,REC=J) LINE(1:2*NPX1)
                DO I=1,NPX1
                  K0=ICHAR(LINE(2*I-1:2*I-1))
                  K1=ICHAR(LINE(2*I:2*I))
                  DN(I,J)=K0+256*K1
                ENDDO
              ENDDO
            CLOSE(UNIT=10)
          ENDIF
          J=NINT(CPIC(2))
          CPIC(1)=CPIC(1)-NINT(Z1)+NINT(DPX(J))
        ENDIF
        NLN1=NLN
      ENDIF

      IF(PICNM(1:1).EQ.'W') THEN
      write (*,*) "Status Www "
        PICFILE='NONE'
        EX=.TRUE.
        DO K0=0,NLN/HGT(KF)-2
          DO I=1,NPX
          DO J=1,2*HGT(KF)
            DN1(I,J)=0
            DN2(I,J)=0
          ENDDO
          ENDDO
          DO I=1,NPX
          DO J=1,HGT(KF)
            DN1(I,J)=DN(I,J+K0*HGT(KF))
            DN2(I,J+HGT(KF))=DN(I,J+(K0+1)*HGT(KF))
          ENDDO
          ENDDO
          K=K0+1
          call correlator(dn1,dn2,npx,2*HGT(KF), x(k),y(k),found(k))
          IF(FOUND(K)) THEN
            WRITE(6,FMT='(I5,2F15.5)') K, X(K), Y(K)
          ELSE
            WRITE(6,FMT='(I5,30X)') K
            X(K)=0
            Y(K)=0
          ENDIF
        ENDDO
        IF(EX) THEN
          DO K=2,NLN/HGT(KF)-1
            X(K)=X(K)+X(K-1)
            Y(K)=Y(K)+Y(K-1)
          ENDDO
          X1= 0
          X2= 0
          DO K=1,NLN/HGT(KF)-1
            x1=min(x1,x(k))
            x2=max(x2,x(k))
          ENDDO
          NPX1=NPX+NINT(X2-X1)
          NLN1=NLN-NINT(Y(NFM-1))
          K=(NFM-1)/2
          CPIC(1)=CPIC(1)+X2-X(K)
          CPIC(2)=K*HGT(KF)+0.5*(HGT(KF)+1)-Y(K)
          DO I=1,NPX1
          DO J=1,NLN1
            DN1(I,J)=0
            DN2(I,J)=0
          ENDDO
          ENDDO
c  ......................................
c    remove hash marks from mono images
c  ......................................
          IF(PICNM(11:11).EQ.'M') THEN
      write (*,*) "Status Mmm "
            DO I=1,NPX
            DO J=1,NLN/HGT(KF)
              DN(I,J*HGT(KF))=0
            ENDDO
            ENDDO
          ENDIF
c  ......................................
          DO I=1,NPX
          DO J=1,HGT(KF)
          IF(DN(I,J).NE.0) THEN
            DN1(I+NINT(X2),J)=DN(I,J)
            DN2(I+NINT(X2),J)=1
          ENDIF
          ENDDO
          ENDDO
          DO K=1,NLN/HGT(KF)-1
            DO I=1,NPX
            DO J=K*HGT(KF)+1,(K+1)*HGT(KF)
            IF(DN(I,J).NE.0) THEN
              DN1(I-NINT(X(K)-X2),J-NINT(Y(K)))=
     .        DN1(I-NINT(X(K)-X2),J-NINT(Y(K)))+DN(I,J)
              DN2(I-NINT(X(K)-X2),J-NINT(Y(K)))=
     .        DN2(I-NINT(X(K)-X2),J-NINT(Y(K)))+1
            ENDIF
            ENDDO
            ENDDO
          ENDDO
          DO I=1,NPX1
          DO J=1,NLN1
            IF(DN2(I,J).NE.0) DN1(I,J)=DN1(I,J)/DN2(I,J)
            DN(I,J)=DN1(I,J)
          ENDDO
          ENDDO
          I=SLEN(PICNM)
          PICFILE='IMAGEFILES/'//PICNM(1:I)//'.DAT'
          IF(BITS.EQ.8) THEN
            OPEN(UNIT=10,FILE=PICFILE,ACCESS='DIRECT',RECL=NPX1,
     .           STATUS='UNKNOWN')
              DO J=1,NLN1
                DO I=1,NPX1
                  LINE(I:I)=CHAR(DN(I,J))
                ENDDO
                WRITE(10,REC=J) LINE(1:NPX1)
              ENDDO
            CLOSE(UNIT=10)
          ELSE
            OPEN(UNIT=10,FILE=PICFILE,ACCESS='DIRECT',RECL=2*NPX1,
     .           STATUS='UNKNOWN')
              DO J=1,NLN1
                DO I=1,NPX1
                  K1=DN(I,J)/256
                  K0=DN(I,J)-256*K1
                  LINE(2*I-1:2*I-1)=CHAR(K0)
                  LINE(2*I:2*I)=CHAR(K1)
                ENDDO
                WRITE(10,REC=J) LINE(1:2*NPX1)
              ENDDO
            CLOSE(UNIT=10)
          ENDIF
        ENDIF
        SEC_PER_LINE=(ET1-ET0)/(NLN1-HGT(KF))
      ENDIF

      DNMIN=1000000
      DNMAX=0
      DO J=1,NLN1
      DO I=1,NPX1
        IF(DN(I,J).NE.0) THEN
          DNMAX=MAX(DNMAX,DN(I,J))
          DNMIN=MIN(DNMIN,DN(I,J))
        ENDIF
      ENDDO
      ENDDO

      OPEN(UNIT=10, FILE='TEMPFILE.gray',RECL=NPX1,
     .     ACCESS='DIRECT', STATUS='UNKNOWN')
        DO J=1,NLN1
          DO I=1,NPX1
            LINE(I:I)=CHAR(0)
            IF(DN(I,J).NE.0) THEN
              LINE(I:I)=CHAR(NINT(1+254.*(DN(I,J)-DNMIN)/(DNMAX-DNMIN)))
            ENDIF
          ENDDO
          WRITE(10,REC=J) LINE(1:NPX1)
        ENDDO
      CLOSE(UNIT=10)
      infile= 'TEMPFILE.gray' 
      tmpfile= 'TEMPFILE1.pgm'  
      call raw2pgm(infile,tmpfile,NPX1,NLN1)
      open(unit=63, file=infile, status='old')
      close(unit=63, status='delete')
      write(6,*) 'gc TEMPFILE1.pgm'

      POL=SEC_PER_LINE*DEG_PER_DAY*RPD()/SPD()

      CALL SPKEZ (10, ET, BFRAME, ABCORR, TARGET, VDUM, LT)
      CALL VEQU(VDUM,USUN)
      CALL VHAT ( USUN, USUN )

      DO I=1,2
      DO J=1,2
        MAT(I,J)=0
      ENDDO
      ENDDO
      DO I=1,2
        W(I)=0
      ENDDO
      DO I=1,NLN1
C        ETL=ET0+(I-0.5)*SEC_PER_LINE
        ETL=ET+(I-CPIC(2))*SEC_PER_LINE
        CALL SPKEZ (TARGET, ETL, BFRAME, ABCORR, OBSRVR, VDUM, LT)
        CALL VEQU(VDUM,V)
        Z1=POL*(I-CPIC(2))
        CALL VROTV(V,POL_AXIS,Z1,V)
        CALL VHAT(V,V)
        MAT(1,1)=MAT(1,1)+V(1)*V(1)
        MAT(1,2)=MAT(1,2)+V(1)*V(2)
        MAT(2,1)=MAT(2,1)+V(2)*V(1)
        MAT(2,2)=MAT(2,2)+V(2)*V(2)
        W(1)=W(1)+V(1)*V(3)
        W(2)=W(2)+V(2)*V(3)
      ENDDO
      CALL INVERTN(2,MAT,IMAT)
      DO I=1,2
        V(I)=0
        DO J=1,2
          V(I)=V(I)+IMAT(I,J)*W(J)
        ENDDO
      ENDDO
      V(3)=-1.d0
      CALL VHAT(V,TRJ_AXIS)

      DO I=1,3
      DO J=1,3
        MAT(I,J)=0
      ENDDO
      ENDDO
      DO I=1,3
        W(I)=0
      ENDDO
      DO I=1,NLN1
        Z2=I-CPIC(2)
C        ETL=ET0+(I-0.5)*SEC_PER_LINE
        ETL=ET+(I-CPIC(2))*SEC_PER_LINE
        CALL SPKEZ (TARGET, ETL, BFRAME, ABCORR, OBSRVR, VDUM, LT)
        CALL VEQU(VDUM,V)
        DO K1=1,3
          DO K2=1,3
            IF((K1.EQ.1).AND.(K2.EQ.1)) THEN
              MAT(K1,K2)=MAT(K1,K2)+1.D0
            ELSE
              MAT(K1,K2)=MAT(K1,K2)+Z2**(K1+K2-2)
            ENDIF
          ENDDO
          IF(K1.EQ.1) THEN
            W(K1)=W(K1)+VNORM(V)
          ELSE
            W(K1)=W(K1)+VNORM(V)*Z2**(K1-1)
          ENDIF
        ENDDO
      ENDDO
      CALL INVERTN(3,MAT,IMAT)
      DO I=1,3
        R0(I)=0
        DO J=1,3
          R0(I)=R0(I)+IMAT(I,J)*W(J)
        ENDDO
      ENDDO

      CALL SPKEZ (TARGET, ET, BFRAME, ABCORR, OBSRVR, VDUM, LT)
      CALL VEQU(VDUM,V0)
      CALL VHAT(V0,U0)

c      Z1=VDOT(TRJ_AXIS,V0)
c      CALL VSCL(Z1,TRJ_AXIS,U0)
c      CALL VSUB(V0,U0,U0)
c      CALL VHAT(U0,U0)

      DO N=1,3
        DO I=1,3
        DO J=1,3
          MAT(I,J)=0
        ENDDO
        ENDDO
        DO I=1,3
          W(I)=0
        ENDDO
        DO I=1,NLN1
          Z2=I-CPIC(2)
          ETL=ET+(I-CPIC(2))*SEC_PER_LINE
          CALL SPKEZ (TARGET, ETL,BFRAME,ABCORR,OBSRVR,VDUM,LT)
          CALL VEQU(VDUM,V)
          Z1=POL*Z2
          CALL VROTV(V,POL_AXIS,Z1,V)
          CALL VHAT(V,V)
          CALL VCRSS(U0,V,V)
          Z1=ASIN(VDOT(TRJ_AXIS,V))
          DO K1=1,3
            DO K2=1,3
              MAT(K1,K2)=MAT(K1,K2)+Z2**(K1+K2)
            ENDDO
            W(K1)=W(K1)+Z1*Z2**K1
          ENDDO
        ENDDO
        CALL INVERTN(3,MAT,IMAT)
        DO I=1,3
          TRJ(I)=0
          DO J=1,3
            TRJ(I)=TRJ(I)+IMAT(I,J)*W(J)
          ENDDO
        ENDDO
        DO I=1,3
          U0(I)=0
        ENDDO
        DO I=1,NLN1
          Z2=I-CPIC(2)
          ETL=ET+(I-CPIC(2))*SEC_PER_LINE
          CALL SPKEZ (TARGET,ETL,BFRAME,ABCORR,OBSRVR,VDUM,LT)
          CALL VEQU(VDUM,V)
          Z1=POL*Z2
          CALL VROTV(V,POL_AXIS,Z1,V)
          Z1=TRJ(1)*Z2+TRJ(2)*Z2**2+TRJ(3)*Z2**3
          CALL VROTV(V,TRJ_AXIS,-Z1,V)
          U0(1)=U0(1)+V(1)
          U0(2)=U0(2)+V(2)
          U0(3)=U0(3)+V(3)
        ENDDO
        CALL VHAT(U0,U0)
      ENDDO
      WRITE(6,FMT='(A12,3D15.6)') 'R0 =        ', (R0(I), I=1,3)
      WRITE(6,FMT='(A12,3D15.6)') 'U0 =        ', (U0(I), I=1,3)
      WRITE(6,FMT='(A12,3D15.6)') 'TRJ_AXIS =  ', (TRJ_AXIS(I), I=1,3)
      WRITE(6,FMT='(A12,3D15.6)') 'TRJ =       ', (TRJ(I), I=1,3)

      Z1=0
      DO I=1,NLN1
        ETL=ET+(I-CPIC(2))*SEC_PER_LINE
        CALL SPKEZ (TARGET, ETL, BFRAME, ABCORR, OBSRVR, VDUM, LT)
        CALL VEQU(VDUM,V)
        Z2=I
        CALL SCOBJ(R0,U0,TRJ_AXIS,TRJ,POL_AXIS,POL,CPIC(2),Z2,W)
        CALL VSUB(V,W,W)
        Z1=Z1+VDOT(W,W)
      ENDDO
      DSCOBJ=SQRT(Z1/NLN1)
      WRITE(6,FMT='(A12,D15.6)') 'DSCOBJ =    ', DSCOBJ

      KMAX=3
40    CONTINUE

      IF(PICNM(1:1).EQ.'W') THEN
        IF(SUF.EQ.'1') Z2=ZW(1) 
        IF(SUF.EQ.'2') Z2=ZW(2)
        IF(SUF.EQ.'3') Z2=ZW(3)
        IF(SUF.EQ.'4') Z2=ZW(4)
        IF(SUF.EQ.'5') Z2=ZW(5)
        IF(SUF.EQ.'6') Z2=ZW(6)
        IF(SUF.EQ.'7') Z2=ZW(7)
        ROTP=-ASIN(Z2/SQRT(Z2**2+MMFL**2))  
        MMFL=SQRT(Z2**2+MMFL**2)              !!!
      ENDIF

      CALL PXFORM (INST, BFRAME, ET, MA)
      DO I=1,3
        CY(I)=-MA(I,1)
        CX(I)= MA(I,2)
        CZ(I)= MA(I,3)
      ENDDO
      IF(PICNM(1:1).EQ.'W') THEN   ! check order
        CALL VROTV(CY,CX,ROTP,CY)        
        CALL VROTV(CZ,CX,ROTP,CZ)        
      ENDIF
      IF(FLIP) THEN
        CALL VMINUS(CX,CX)
        CALL VMINUS(CY,CY)
      ENDIF
c      IF(PICNM(1:1).EQ.'W') THEN   ! check order
c        CALL VROTV(CY,CX,ROTP,CY)        
c        CALL VROTV(CZ,CX,ROTP,CZ)        
c      ENDIF
      WRITE(6,FMT='(A12,2F15.4)') 'CPIC =      ', (CPIC(I), I=1,2)
      WRITE(6,FMT='(A12,3D15.6)') 'CX =        ', (CX(I), I=1,3)
      WRITE(6,FMT='(A12,3D15.6)') 'CY =        ', (CY(I), I=1,3)
      WRITE(6,FMT='(A12,3D15.6)') 'CZ =        ', (CZ(I), I=1,3)

      DO N=1,3        
        DO I=1,KMAX+1
        DO J=1,KMAX+1
          MAT(I,J)=0
        ENDDO
        ENDDO
        DO J=1,3
        DO I=1,KMAX+1
          WMAT(I,J)=0.D0
        ENDDO
        ENDDO
        DO I=1,NLN1
          ETL=ET+(I-CPIC(2))*SEC_PER_LINE
          CALL PXFORM (INST, BFRAME, ETL, MA)
          DO J=1,3
            VY(J)=-MA(J,1)
            VX(J)= MA(J,2)
            VZ(J)= MA(J,3)
          ENDDO
          IF(PICNM(1:1).EQ.'W') THEN   ! check order
            CALL VROTV(VY,VX,ROTP,VY)        
            CALL VROTV(VZ,VX,ROTP,VZ)        
          ENDIF
          IF(FLIP) THEN
            CALL VMINUS(VX,VX)
            CALL VMINUS(VY,VY)
          ENDIF
c          IF(PICNM(1:1).EQ.'W') THEN   ! check order
c            CALL VROTV(VY,VX,ROTP,VY)        
c            CALL VROTV(VZ,VX,ROTP,VZ)        
c          ENDIF
          Z2=I-CPIC(2)
          Z1=POL*Z2
          CALL VROTV(VX,POL_AXIS,Z1,VX)
          CALL VROTV(VY,POL_AXIS,Z1,VY)
          CALL VROTV(VZ,POL_AXIS,Z1,VZ)
          Z1=TRJ(1)*Z2+TRJ(2)*Z2**2+TRJ(3)*Z2**3
          CALL VROTV(VX,TRJ_AXIS,-Z1,VX)
          CALL VROTV(VY,TRJ_AXIS,-Z1,VY)
          CALL VROTV(VZ,TRJ_AXIS,-Z1,VZ)
          MA(1,1)=VDOT(VX,CX)
          MA(1,2)=VDOT(VX,CY)
          MA(1,3)=VDOT(VX,CZ)
          MA(2,1)=VDOT(VY,CX)
          MA(2,2)=VDOT(VY,CY)
          MA(2,3)=VDOT(VY,CZ)
          MA(3,1)=VDOT(VZ,CX)
          MA(3,2)=VDOT(VZ,CY)
          MA(3,3)=VDOT(VZ,CZ)
          CALL MAT2ROT(MA,V)
          DO K1=1,KMAX+1
            DO K2=1,KMAX+1
              IF((K1.EQ.1).AND.(K2.EQ.1)) THEN
                MAT(K1,K2)=MAT(K1,K2)+1.D0
              ELSE
                MAT(K1,K2)=MAT(K1,K2)+Z2**(K1+K2-2)
              ENDIF
            ENDDO
            DO J=1,3
              IF(K1.EQ.1) THEN
                WMAT(K1,J)=WMAT(K1,J)+V(J)
              ELSE
                WMAT(K1,J)=WMAT(K1,J)+V(J)*Z2**(K1-1)
              ENDIF
            ENDDO
          ENDDO
        ENDDO
        CALL INVERTN(KMAX+1,MAT,IMAT)
        DO K=1,3
        DO I=1,KMAX+1
          CMAT(I,K)=0.D0
          DO J=1,KMAX+1
            CMAT(I,K)=CMAT(I,K)+IMAT(I,J)*WMAT(J,K)
          ENDDO
        ENDDO
        ENDDO
        DO J=1,3
          V(J)=CMAT(1,J)
        ENDDO
        CALL ROT2MAT(V,MB)
        DO J=1,3
          MA(1,J)=CX(J)
          MA(2,J)=CY(J)
          MA(3,J)=CZ(J)
        ENDDO
        CALL MXM(MB,MA,MA)
        DO J=1,3
          CX(J)=MA(1,J)
          CY(J)=MA(2,J)
          CZ(J)=MA(3,J)
        ENDDO
      ENDDO
      DO I=1,KMAX
      DO J=1,3
        CAX(I,J)=CMAT(I+1,J)
      ENDDO
      ENDDO
      WRITE(6,FMT='(A12,3D15.6)') 'CX =        ', (CX(I), I=1,3)
      WRITE(6,FMT='(A12,3D15.6)') 'CY =        ', (CY(I), I=1,3)
      WRITE(6,FMT='(A12,3D15.6)') 'CZ =        ', (CZ(I), I=1,3)
      WRITE(6,FMT='(A12,3D15.6)') 'SZ =        ', (USUN(I), I=1,3)
      DO K=1,KMAX
        WRITE(6,FMT='(A12,3D15.6)') 'CAX =       ', (CAX(K,I), I=1,3)
      ENDDO

      V0(1)=0
      V0(2)=0
      V0(3)=0
      V(1)=1
      V(2)=1
      V(3)=1
      DO I=1,NLN1
        ETL=ET+(I-CPIC(2))*SEC_PER_LINE
        CALL PXFORM (INST, BFRAME, ETL, MA)
        DO J=1,3
          VY(J)=-MA(J,1)
          VX(J)= MA(J,2)
          VZ(J)= MA(J,3)
        ENDDO
        IF(PICNM(1:1).EQ.'W') THEN   ! check order
          CALL VROTV(VY,VX,ROTP,VY)        
          CALL VROTV(VZ,VX,ROTP,VZ)        
        ENDIF
        IF(FLIP) THEN
          CALL VMINUS(VX,VX)
          CALL VMINUS(VY,VY)
        ENDIF
c        IF(PICNM(1:1).EQ.'W') THEN   ! check order
c          CALL VROTV(VY,VX,ROTP,VY)        
c          CALL VROTV(VZ,VX,ROTP,VZ)        
c        ENDIF
        Z1=I
        CALL POINTING(CX,CY,CZ,CAX,KMAX,
     .                TRJ_AXIS,TRJ,POL_AXIS,POL,CPIC(2),Z1,WX,WY,WZ)
        V0(1)=V0(1)+VDOT(VX,WX)
        V0(2)=V0(2)+VDOT(VY,WY)
        V0(3)=V0(3)+VDOT(VZ,WZ)
        V(1)=MIN(V(1),VDOT(VX,WX))
        V(2)=MIN(V(2),VDOT(VY,WY))
        V(3)=MIN(V(3),VDOT(VZ,WZ))
      ENDDO
      DO I=1,3
        DPTG(I)=ACOS(V0(I)/NLN1)
      ENDDO
      WRITE(6,FMT='(A12,3D15.6)') 'DPTG =      ', 
     .                            (DPTG(I), I=1,3)
      WRITE(6,FMT='(A12,3D15.6)') 'MXDPTG =    ', 
     .                            (ACOS(V(I)), I=1,3)
      Z1=1/(K11*MMFL)
      WRITE(6,FMT='(A12,D15.6)') 'PXSIZE =    ', Z1
      WRITE(6,FMT='(A12,3D15.6)') 'U0CAM =     ', VDOT(U0,CX),
     .                             VDOT(U0,CY), VDOT(U0,CZ)
      WRITE(6,*) 
      WRITE(6,*) 'KMAX = ', KMAX
      WRITE(6,*) 'Change KMAX < 20? (y/n)'
      READ(5,FMT='(A1)') ANS
      IF(ANS.EQ.'y') THEN
        WRITE(6,*) 'Input new KMAX'
        READ(5,*) KMAX
        GO TO 40
      ENDIF

      Z2=(1.D0+NLN1)/2
      CALL SCOBJ(R0,U0,TRJ_AXIS,TRJ,POL_AXIS,POL,CPIC(2),Z2,W)
      CALL VMINUS(W,U)
      CALL U2V(U,V)
      R00=SQRT(VDOT(V,V))
      write(6,*) r00          !!!!

      DO I=1,8
      DO J=1,8
        MAT(I,J)=0
      ENDDO
      ENDDO
      DO I=1,8
        WW(I)=0
        XX(I)=0
      ENDDO
c      XX(1)=CPX(PICNM,NPX1,NLN1,CPIC,NINT(CPIC(2)))
      XX(1)=CPIC(1)
      XX(2)=K11
      XX(5)=CPIC(2)
      MAT(1,1)=1.D20
      MAT(2,2)=1.D20
      MAT(5,5)=1.D20
      CALL VSCL(R0(1),U0,V0)
      DO J=1,NLN1,20
        X2=CPX(PICNM,NPX1,NLN1,CPIC,J)              
        Z2=J
        CALL SCOBJ(R0,U0,TRJ_AXIS,TRJ,POL_AXIS,POL,CPIC(2),Z2,W)
        DO I=1,NPX1,20
          Z1=I
          Z2=J
          IF(ONPX(PICNM,NPX1,NLN1,I,J)) THEN
            CALL PL2U(CX,CY,CZ,CAX,KMAX,TRJ_AXIS,TRJ,
     .              POL_AXIS,POL,MMFL,K11,X2,CPIC(2),Z1,Z2,U)
            Z1=VDOT(U,W)
            Z2=Z1-SQRT(Z1**2+R00**2-VDOT(W,W))
            CALL VSCL(Z2,U,V)
            CALL VSUB(V,W,V)
            Z1=VDOT(CX,V)+VDOT(CX,V0)
            Z2=VDOT(CY,V)+VDOT(CY,V0)
            Z3=VDOT(CZ,V)+VDOT(CZ,V0)
            X1=MMFL*Z1/Z3
            Y1=MMFL*Z2/Z3
            Z1=Y1/MMFL
            Z2=VDOT(U0,CY)
            Z3=SQRT(1-Z2**2)
            Y1=Y1/(Z3+Z1*Z2)
            PAR(1)=1
            PAR(2)=X1
            PAR(3)=Y1
            PAR(4)=X1*Y1
            DO K1=1,4
            DO K2=1,4
              MAT(K1,K2)=MAT(K1,K2)+PAR(K1)*PAR(K2)            
              MAT(K1+4,K2+4)=MAT(K1+4,K2+4)+PAR(K1)*PAR(K2)            
            ENDDO
            ENDDO
            Z1=0
            Z2=0
            DO K1=1,4
              Z1=Z1+PAR(K1)*XX(K1)
              Z2=Z2+PAR(K1)*XX(K1+4)
            ENDDO
            DO K1=1,4
              WW(K1)=WW(K1)+PAR(K1)*(I-Z1)
              WW(K1+4)=WW(K1+4)+PAR(K1)*(J-Z2)
            ENDDO
          ENDIF
        ENDDO
      ENDDO
      WRITE(6,*) K11
      CALL INVERTN(8,MAT,IMAT)
      DO I=1,8
        XX(I)=0
      ENDDO
c      XX(1)=CPX(PICNM,NPX1,NLN1,CPIC,NINT(CPIC(2)))
      XX(1)=CPIC(1)
      XX(2)=K11
      XX(5)=CPIC(2)
      DO I=1,8
      DO J=1,8
        XX(I)=XX(I)+IMAT(I,J)*WW(J)
      ENDDO
      ENDDO
      WRITE(6,FMT='(4F15.4)') (XX(I), I=1,4)
      WRITE(6,FMT='(4F15.4)') (XX(I+4), I=1,4)
      X3=0
      Y3=0
      K=0
      DO J=1,NLN1,20
        X2=CPX(PICNM,NPX1,NLN1,CPIC,J)              
        Z2=J
        CALL SCOBJ(R0,U0,TRJ_AXIS,TRJ,POL_AXIS,POL,CPIC(2),Z2,W)
        DO I=1,NPX1,20
          Z1=I
          Z2=J
          IF(ONPX(PICNM,NPX1,NLN1,I,J)) THEN
            CALL PL2U(CX,CY,CZ,CAX,KMAX,TRJ_AXIS,TRJ,
     .              POL_AXIS,POL,MMFL,K11,X2,CPIC(2),Z1,Z2,U)
            Z1=VDOT(U,W)
            Z2=Z1-SQRT(Z1**2+R00**2-VDOT(W,W))
            CALL VSCL(Z2,U,V)
            CALL VSUB(V,W,V)
            Z1=VDOT(CX,V)+VDOT(CX,V0)
            Z2=VDOT(CY,V)+VDOT(CY,V0)
            Z3=VDOT(CZ,V)+VDOT(CZ,V0)
            X1=MMFL*Z1/Z3
            Y1=MMFL*Z2/Z3
            Z1=Y1/MMFL
            Z2=VDOT(U0,CY)
            Z3=SQRT(1-Z2**2)
            Y1=Y1/(Z3+Z1*Z2)
            PAR(1)=1
            PAR(2)=X1
            PAR(3)=Y1
            PAR(4)=X1*Y1
            Z1=0
            Z2=0
            DO K1=1,4
              Z1=Z1+PAR(K1)*XX(K1)
              Z2=Z2+PAR(K1)*XX(K1+4)
            ENDDO
            X3=X3+(I-Z1)**2
            Y3=Y3+(J-Z2)**2
            K=K+1
          ENDIF
        ENDDO
      ENDDO
      WRITE(6,*) SQRT(X3/K), SQRT(Y3/K), SQRT((X3+Y3)/K)

      OPEN(UNIT=25, FILE='RESIDUALS.TXT', STATUS='UNKNOWN')
      DO J=1,NLN1,100
        X2=CPX(PICNM,NPX1,NLN1,CPIC,J)              
        Z2=J
        CALL SCOBJ(R0,U0,TRJ_AXIS,TRJ,POL_AXIS,POL,CPIC(2),Z2,W)
        DO I=1,NPX1,100
          Z1=I
          Z2=J
          IF(ONPX(PICNM,NPX1,NLN1,I,J)) THEN
            CALL PL2U(CX,CY,CZ,CAX,KMAX,TRJ_AXIS,TRJ,
     .              POL_AXIS,POL,MMFL,K11,X2,CPIC(2),Z1,Z2,U)
            Z1=VDOT(U,W)
            Z2=Z1-SQRT(Z1**2+R00**2-VDOT(W,W))
            CALL VSCL(Z2,U,V)
            CALL VSUB(V,W,V)
            Z1=VDOT(CX,V)+VDOT(CX,V0)
            Z2=VDOT(CY,V)+VDOT(CY,V0)
            Z3=VDOT(CZ,V)+VDOT(CZ,V0)
            X1=MMFL*Z1/Z3
            Y1=MMFL*Z2/Z3
            Z1=Y1/MMFL
            Z2=VDOT(U0,CY)
            Z3=SQRT(1-Z2**2)
            Y1=Y1/(Z3+Z1*Z2)
            PAR(1)=1
            PAR(2)=X1
            PAR(3)=Y1
            PAR(4)=X1*Y1
            Z1=0
            Z2=0
            DO K1=1,4
              Z1=Z1+PAR(K1)*XX(K1)
              Z2=Z2+PAR(K1)*XX(K1+4)
            ENDDO
            WRITE(25,FMT='(2I10,2F10.2)') I, J, I-Z1, J-Z2
          ENDIF
        ENDDO
      ENDDO
      CLOSE(UNIT=25)

      WRITE(6,*) 'Keep image? (y/n)'
      READ(5,FMT='(A1)') ANS
      IF(ANS.NE.'y') THEN
        WRITE(6,*) 'This deletes existing SUMFILE.  Continue? (y/n)'
        READ(5,FMT='(A1)') ANS
        IF(ANS.EQ.'y') THEN
          OPEN(UNIT=10, FILE=PICFILE, STATUS='OLD')
          CLOSE(UNIT=10, STATUS='DELETE')
          I=SLEN(PICNM)
          PICFILE='SUMFILES/'//PICNM(1:I)//'.SUM'
          OPEN(UNIT=10, FILE=PICFILE, STATUS='UNKNOWN')
          CLOSE(UNIT=10, STATUS='DELETE')
          PICFILE='NOMINALS/'//PICNM(1:I)//'.NOM'
          OPEN(UNIT=10, FILE=PICFILE, STATUS='UNKNOWN')
          CLOSE(UNIT=10, STATUS='DELETE')
        ENDIF
      ELSE
        ANS='a'
        I=SLEN(PICNM)
        PICFILE='SUMFILES/'//PICNM(1:I)//'.SUM'
        INQUIRE(FILE=PICFILE, EXIST=EX)
        IF(EX) THEN
          WRITE(6,*) 'a. Overwrite SUMFILE'
          WRITE(6,*) 'b. Change geometry only'
          WRITE(6,*) 'c. No Changes'
          READ(5,FMT='(A1)') ANS
        ENDIF
        IF(ANS.EQ.'a') THEN
          OPEN(UNIT=20,FILE=PICFILE,STATUS='UNKNOWN')
            WRITE(20,FMT='(A12,1X,A1,66X)') PICNM, 'L'
            WRITE(20,FMT='(A24,56X)') UTC
            WRITE(20,FMT='(4I6,36X,A20)') NPX1, NLN1, T1, T2,
     .                                  '  NPX, NLN, THRSH  '
            WRITE(20,FMT='(3D20.10,A20)') MMFL, XX(1), XX(5),
     .                                  '  MMFL, CTR        '
            CALL VSCL(R0(1),U0,V0)
            WRITE(20,FMT='(3D20.10,A20)') (V0(I), I=1,3), 
     .                                  '  SCOBJ            '
            WRITE(20,FMT='(3D20.10,A20)') (CX(I), I=1,3),
     .                                  '  CX               ' 
            WRITE(20,FMT='(3D20.10,A20)') (CY(I), I=1,3),
     .                                  '  CY               ' 
            WRITE(20,FMT='(3D20.10,A20)') (CZ(I), I=1,3),
     .                                  '  CZ               ' 
            WRITE(20,FMT='(3D20.10,A20)') (USUN(I), I=1,3),
     .                                  '  SZ               ' 
            WRITE(20,FMT='(6F10.4,A20)') (XX(I), I=2,4), (XX(I), I=6,8), 
     .                                  '  K-MATRIX         '
            WRITE(20,FMT='(4D15.5,A20)') (D(I), I=1,4),
     .                                  '  DISTORTION       '
            W(1)=SQRT(DSCOBJ**2+SIG_VSO(1)**2)
            W(2)=SQRT(DSCOBJ**2+SIG_VSO(2)**2)
            W(3)=SQRT(DSCOBJ**2+SIG_VSO(3)**2)
            WRITE(20,FMT='(3D20.10,A20)') (W(I), I=1,3),
     .                                  '  SIGMA_VSO        ' 
            W(1)=SQRT(DPTG(1)**2+SIG_PTG(1)**2)
            W(2)=SQRT(DPTG(2)**2+SIG_PTG(2)**2)
            W(3)=SQRT(DPTG(3)**2+SIG_PTG(3)**2)
            WRITE(20,FMT='(3D20.10,A20)') (W(I), I=1,3),
     .                                  '  SIGMA_PTG        ' 
            WRITE(20,FMT='(A14,66X)')  'LINE SCAN DATA'
            WRITE(20,FMT='(3D20.10,A20)') (U0(I), I=1,3), 
     .                                '  U0               '
            WRITE(20,FMT='(3D20.10,A20)') (R0(I), I=1,3),
     .                                  '  R0               ' 
            WRITE(20,FMT='(3D20.10,A20)') (POL_AXIS(I), I=1,3),
     .                                  '  POL_AXIS         '
            WRITE(20,FMT='(D20.10,40X,A20)') POL,
     .                                  '  POL              '
            WRITE(20,FMT='(3D20.10,A20)') (TRJ_AXIS(I), I=1,3),
     .                                  '  TRJ_AXIS         '
            WRITE(20,FMT='(3D20.10,A20)') (TRJ(I), I=1,3),
     .                                  '  TRJ              '
            WRITE(20,FMT='(I5,55X,A20)') KMAX,
     .                                  '  KMAX             '
            DO K=1,KMAX
              WRITE(20,FMT='(3D20.10,A20)') (CAX(K,I), I=1,3),
     .                                  '  CAX              ' 
            ENDDO
            WRITE(20,FMT='(A9,71X)')  'LANDMARKS'
            WRITE(20,FMT='(A9,71X)')  'LIMB FITS'
            WRITE(20,FMT='(A8,72X)')  'END FILE'
          CLOSE(UNIT=20)
        ENDIF
        IF(ANS.EQ.'b') THEN
          OPEN(UNIT=20,FILE=PICFILE,STATUS='OLD')
            I=0
50          CONTINUE
            READ(20,FMT='(A80)') ZLINE0
            IF(ZLINE0(1:9).NE.'LANDMARKS') GO TO 50
60          CONTINUE
            I=I+1
            READ(20,FMT='(A80)') ZLINE(I)
            IF(ZLINE(I)(1:8).NE.'END FILE') GO TO 60
          CLOSE(UNIT=20)
          IMAX=I
          OPEN(UNIT=20,FILE=PICFILE,STATUS='UNKNOWN')
            WRITE(20,FMT='(A12,1X,A1,66X)') PICNM, 'L'
            WRITE(20,FMT='(A24,56X)') UTC
            WRITE(20,FMT='(4I6,36X,A20)') NPX1, NLN1, T1, T2,
     .                                  '  NPX, NLN, THRSH  '
            WRITE(20,FMT='(3D20.10,A20)') MMFL, XX(1), XX(5),
     .                                  '  MMFL, CTR        '
            CALL VSCL(R0(1),U0,V0)
            WRITE(20,FMT='(3D20.10,A20)') (V0(I), I=1,3), 
     .                                  '  SCOBJ            '
            WRITE(20,FMT='(3D20.10,A20)') (CX(I), I=1,3),
     .                                  '  CX               ' 
            WRITE(20,FMT='(3D20.10,A20)') (CY(I), I=1,3),
     .                                  '  CY               ' 
            WRITE(20,FMT='(3D20.10,A20)') (CZ(I), I=1,3),
     .                                  '  CZ               ' 
            WRITE(20,FMT='(3D20.10,A20)') (USUN(I), I=1,3),
     .                                  '  SZ               ' 
            WRITE(20,FMT='(6F10.4,A20)') (XX(I), I=2,4), (XX(I), I=6,8), 
     .                                  '  K-MATRIX         '
            WRITE(20,FMT='(4D15.5,A20)') (D(I), I=1,4),
     .                                  '  DISTORTION       '
            W(1)=SQRT(DSCOBJ**2+SIG_VSO(1)**2)
            W(2)=SQRT(DSCOBJ**2+SIG_VSO(2)**2)
            W(3)=SQRT(DSCOBJ**2+SIG_VSO(3)**2)
            WRITE(20,FMT='(3D20.10,A20)') (W(I), I=1,3),
     .                                  '  SIGMA_VSO        ' 
            W(1)=SQRT(DPTG(1)**2+SIG_PTG(1)**2)
            W(2)=SQRT(DPTG(2)**2+SIG_PTG(2)**2)
            W(3)=SQRT(DPTG(3)**2+SIG_PTG(3)**2)
            WRITE(20,FMT='(3D20.10,A20)') (W(I), I=1,3),
     .                                  '  SIGMA_PTG        ' 
            WRITE(20,FMT='(A14,66X)')  'LINE SCAN DATA'
            WRITE(20,FMT='(3D20.10,A20)') (U0(I), I=1,3), 
     .                                '  U0               '
            WRITE(20,FMT='(3D20.10,A20)') (R0(I), I=1,3),
     .                                  '  R0               ' 
            WRITE(20,FMT='(3D20.10,A20)') (POL_AXIS(I), I=1,3),
     .                                  '  POL_AXIS         '
            WRITE(20,FMT='(D20.10,40X,A20)') POL,
     .                                  '  POL              '
            WRITE(20,FMT='(3D20.10,A20)') (TRJ_AXIS(I), I=1,3),
     .                                  '  TRJ_AXIS         '
            WRITE(20,FMT='(3D20.10,A20)') (TRJ(I), I=1,3),
     .                                  '  TRJ              '
            WRITE(20,FMT='(I5,55X,A20)') KMAX,
     .                                  '  KMAX             '
            DO K=1,KMAX
              WRITE(20,FMT='(3D20.10,A20)') (CAX(K,I), I=1,3),
     .                                  '  CAX              ' 
            ENDDO
            WRITE(20,FMT='(A9,71X)')  'LANDMARKS'
            DO I=1,IMAX
              WRITE(20,FMT='(A80)')  ZLINE(I)
            ENDDO
          CLOSE(UNIT=20)
        ENDIF
        IF((ANS.NE.'a').AND.(ANS.NE.'b')) GO TO 220

        I=SLEN(PICNM)
        PICFILE='NOMINALS/'//PICNM(1:I)//'.NOM'
        OPEN(UNIT=20,FILE=PICFILE,STATUS='UNKNOWN')
          WRITE(20,FMT='(A12,68X)') PICNM
          WRITE(20,FMT='(A24,56X)') UTC
          WRITE(20,FMT='(3D20.10,A20)') (V0(I), I=1,3), 
     .                                '  SCOBJ            '
          W(1)=SQRT(DSCOBJ**2+SIG_VSO(1)**2)
          W(2)=SQRT(DSCOBJ**2+SIG_VSO(2)**2)
          W(3)=SQRT(DSCOBJ**2+SIG_VSO(3)**2)
          WRITE(20,FMT='(3D20.10,A20)') (W(I), I=1,3),
     .                                '  SIGMA_VSO        ' 
          WRITE(20,FMT='(3D20.10,A20)') (CX(I), I=1,3),
     .                                '  CX               ' 
          WRITE(20,FMT='(3D20.10,A20)') (CY(I), I=1,3),
     .                                '  CY               ' 
          WRITE(20,FMT='(3D20.10,A20)') (CZ(I), I=1,3),
     .                                '  CZ               ' 
          W(1)=SQRT(DPTG(1)**2+SIG_PTG(1)**2)
          W(2)=SQRT(DPTG(2)**2+SIG_PTG(2)**2)
          W(3)=SQRT(DPTG(3)**2+SIG_PTG(3)**2)
          WRITE(20,FMT='(3D20.10,A20)') (W(I), I=1,3),
     .                                '  SIGMA_PTG        ' 
          WRITE(20,FMT='(A8,72X)')  'END FILE'
        CLOSE(UNIT=20)

        INQUIRE(FILE='PICTLIST.TXT', EXIST=EX)
        IF(.NOT.EX) THEN
          OPEN(UNIT=20,FILE='PICTLIST.TXT',STATUS='NEW')
            WRITE(20,FMT='(A3,10X)') 'END'        
          CLOSE(UNIT=20)
        ENDIF

        I=0
        OPEN(UNIT=20,FILE='PICTLIST.TXT',STATUS='OLD')
120       CONTINUE
          I=I+1
          READ(20,FMT='(A13)') XNAME(I)
          IF(XNAME(I)(2:13).EQ.PICNM) THEN
            CLOSE(UNIT=20)
            GO TO 220
          ENDIF
          IF(XNAME(I)(1:3).NE.'END') GO TO 120
        CLOSE(UNIT=20)
        OPEN(UNIT=20,FILE='PICTLIST.TXT',STATUS='OLD')
          DO J=1,I-1  
            WRITE(20,FMT='(A13)') XNAME(J)
          ENDDO
          WRITE(20,FMT='(1X,A12)') PICNM        
          WRITE(20,FMT='(A3,10X)') 'END'        
        CLOSE(UNIT=20)

220     CONTINUE
      ENDIF

      OPEN(UNIT=10,FILE='TEMPFILE0.pgm',STATUS='UNKNOWN')
      CLOSE(UNIT=10,STATUS='DELETE')
      OPEN(UNIT=10,FILE='TEMPFILE1.pgm',STATUS='UNKNOWN')
      CLOSE(UNIT=10,STATUS='DELETE')

      WRITE(6,*) 'Add more images? (y/n)'
      READ(5,FMT='(A1)') ANS
      IF(ANS.EQ.'y') go to 30

      STOP
      END
      
c  ..............................................
      subroutine raw2pgm(infile,outfile,npx,nln)
c  ..............................................

      implicit none
      
      integer*4         npx, nln, k
      character*20      header
      character*72      infile, outfile
      character*99999   cline0, cline1
      
      do k=1,20
        header(k:k)=' '
      enddo
      header(1:2)='P5'
      write(header(5:9),fmt='(i5)') npx
      write(header(11:15),fmt='(i5)') nln
      header(17:20)='255 '

      open(unit=20, file=outfile,status='unknown')
      close(unit=20, status='delete')
      
      open(unit=20, file=outfile, recl=npx, access='direct', 
     .     status='unknown')
      open(unit=10, file=infile, recl=npx, access='direct', 
     .     status='old')
        do k=1,nln
          read(10,rec=k) cline0(1:npx)
          cline1(1:20)=header
          cline1(21:npx)=cline0(1:npx-20)
          write(20,rec=k) cline1(1:npx)          
          header=cline0(npx-19:npx)
        enddo
        cline1(1:20)=header
        do k=21,npx
          cline1(k:k)=' '
        enddo
        write(20,rec=nln+1) cline1(1:npx)          
      close(unit=10)
      close(unit=20)
      
      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   ..................................................
      FUNCTION NLF(STRING,I)
C   ..................................................

      IMPLICIT NONE

      CHARACTER*(*)       STRING
      INTEGER             I, J, NLF

      J=I
10    J=J+1
      IF(ICHAR(STRING(J:J)).NE.10) GO TO 10
      NLF=J

      RETURN
      END

C   ..................................................
      FUNCTION NEQ(STRING,I)
C   ..................................................

      IMPLICIT NONE

      CHARACTER*(*)       STRING
      INTEGER             I, J, NEQ

      J=I
10    J=J+1
      IF(STRING(J:J).NE.'=') GO TO 10
      NEQ=J

      RETURN
      END

C   ..................................................
      FUNCTION CPX(PICNM,NPX,NLN,CPIC,L)
C   ..................................................

      IMPLICIT NONE

      DOUBLE PRECISION    CPX
      DOUBLE PRECISION    CPIC(2)
      INTEGER             NPX
      INTEGER             NLN
      INTEGER             L
      INTEGER             I
      INTEGER             J
      INTEGER             BYTES
      CHARACTER*12        PICNM
      CHARACTER*72        PICFILE
      CHARACTER*50000     LINE
      LOGICAL             EX

      CALL FIND_PICFILE(PICNM,PICFILE,EX)

      bytes=2
      open(unit=10, file=picfile, access='direct',
     .     recl=2*npx, status='old')
        read(10,rec=nln, err=100) line(1:2*npx)
        go to 101
100     bytes=1
101     continue
      close(unit=10)

      IF(BYTES.EQ.1) THEN
        OPEN(UNIT=22, FILE=PICFILE, ACCESS='DIRECT',
     .       RECL=NPX, STATUS='OLD')
          READ(22,REC=NINT(CPIC(2))) LINE(1:NPX)
          I=0
10        CONTINUE
          I=I+1
          IF(LINE(I:I).EQ.CHAR(0)) GO TO 10
          READ(22,REC=L) LINE(1:NPX)
          J=0
20        CONTINUE
          J=J+1
          IF(LINE(J:J).EQ.CHAR(0)) GO TO 20
        CLOSE(UNIT=22)
      ELSE
        OPEN(UNIT=22, FILE=PICFILE, ACCESS='DIRECT',
     .       RECL=2*NPX, STATUS='OLD')
          READ(22,REC=NINT(CPIC(2))) LINE(1:2*NPX)
          I=0
30        CONTINUE
          I=I+1
          IF((LINE(2*I-1:2*I-1).EQ.CHAR(0)).AND.
     .       (LINE(2*I:2*I).EQ.CHAR(0))) GO TO 30
          READ(22,REC=L) LINE(1:2*NPX)
          J=0
40        CONTINUE
          J=J+1
          IF((LINE(2*J-1:2*J-1).EQ.CHAR(0)).AND.
     .       (LINE(2*J:2*J).EQ.CHAR(0))) GO TO 40
        CLOSE(UNIT=22)
      ENDIF

      CPX=CPIC(1)+J-I

      RETURN
      END

C   ..................................................
      FUNCTION ONPX(PICNM,NPX,NLN,P,L)
C   ..................................................

      IMPLICIT NONE

      INTEGER             NPX
      INTEGER             NLN
      INTEGER             P
      INTEGER             L
      INTEGER             BYTES
      CHARACTER*12        PICNM
      CHARACTER*72        PICFILE
      CHARACTER*50000     LINE
      LOGICAL             ONPX
      LOGICAL             EX

      ONPX=.FALSE.
      IF(P.LT.1) RETURN
      IF(P.GT.NPX) RETURN
      IF(L.LT.1) RETURN
      IF(L.GT.NLN) RETURN
      CALL FIND_PICFILE(PICNM,PICFILE,EX)
      IF(.NOT.EX) RETURN

      bytes=2
      open(unit=10, file=picfile, access='direct',
     .     recl=2*npx, status='old')
        read(10,rec=nln, err=100) line(1:2*npx)
        go to 101
100     bytes=1
101     continue
      close(unit=10)

      if(bytes.eq.1) then
        OPEN(UNIT=22, FILE=PICFILE, ACCESS='DIRECT',
     .       RECL=NPX, STATUS='OLD')
          READ(22,REC=L) LINE(1:NPX)
        CLOSE(UNIT=22)
        IF(LINE(P:P).NE.CHAR(0)) ONPX=.TRUE.
      else
        OPEN(UNIT=22, FILE=PICFILE, ACCESS='DIRECT',
     .       RECL=2*NPX, STATUS='OLD')
          READ(22,REC=L) LINE(1:2*NPX)
        CLOSE(UNIT=22)
        IF((LINE(2*P-1:2*P-1).NE.CHAR(0)).OR.
     .     (LINE(2*P:2*P).NE.CHAR(0))) ONPX=.TRUE.
      endif

      RETURN
      END

c  ...............................................................
      subroutine correlator(dn1,dn2,npx,nln, x0,y0,found)
c  ...............................................................

      implicit none

      integer*4       i, j, i0, j0, i1, j1, i2, k0, imx, jmx
      INTEGER         DN1(2000,20000)
      INTEGER         DN2(2000,20000)
      INTEGER          NPX
      INTEGER          NLN
      
      real*8          z1, z2, z3, z4, cor(-10:10,-10:10), x0, y0, z0 
      
      logical         use(-10:10,-10:10)
      logical         found

        x0=0        
        y0=0        
        z0=-1
        imx=10
        jmx=10
        found=.false.

        do i=-8,8
        do j=-8,8
          z1=0
          z2=0
          z3=0
          z4=0
          i2=0
          cor(i,j)=0
          use(i,j)=.false.
          do i1=1,npx
          do j1=1,nln
            if((i+i1.ge.1).and.(i+i1.le.npx).and.
     .         (j+j1.ge.1).and.(j+j1.le.nln)) then
            if((dn1(i1,j1).gt.0).and.(dn2(i1+i,j1+j).gt.0)) then
              cor(i,j)=cor(i,j)+dn1(i1,j1)*dn2(i1+i,j1+j)
              z1=z1+dn1(i1,j1)
              z2=z2+dn2(i1+i,j1+j)
              z3=z3+dn1(i1,j1)**2
              z4=z4+dn2(i1+i,j1+j)**2
              i2=i2+1
            endif
            endif
          enddo
          enddo
          if(i2.gt.0) then
            cor(i,j)=cor(i,j)/i2
            z1=z1/i2
            z2=z2/i2
            z3=z3/i2
            z4=z4/i2
            z4=(z3-z1**2)*(z4-z2**2)
            if(z4.gt.0) then
              cor(i,j)=(cor(i,j)-z1*z2)/sqrt(z4)
              use(i,j)=.true.
              if(cor(i,j).gt.z0) then
                z0=cor(i,j)
                imx=i
                jmx=j
              endif
            endif
          endif
        enddo
        enddo
        if((abs(imx).le.8).and.(abs(jmx).le.8)) then
          x0=imx
          y0=jmx 
          z1=0
          z2=0
          if(use(imx-1,jmx).and.use(imx+1,jmx)) then
            z1=2*(cor(imx+1,jmx)+cor(imx-1,jmx)-2*cor(imx,jmx))
            if(z1.ne.0) then
              x0=imx+(cor(imx-1,jmx)-cor(imx+1,jmx))/z1
            endif
          endif
          if(use(imx,jmx-1).and.use(imx,jmx+1)) then
            z2=2*(cor(imx,jmx+1)+cor(imx,jmx-1)-2*cor(imx,jmx))
            if(z2.ne.0) then
              y0=jmx+(cor(imx,jmx-1)-cor(imx,jmx+1))/z2
            endif
          endif
          if((z1.ne.0).and.(z2.ne.0)) found=.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 INVERTN(N,M,MINV)
c   ........................................................

      IMPLICIT NONE
      INTEGER*4    J, K, L, N
      REAL*8       M(20,20), MINV(20,20)
      REAL*8       U(20,20), V(20,20), D(20)

      DO K=1,N
      DO L=1,N
        U(K,L)=0.
      ENDDO
      ENDDO

      DO L=1,N
       U(L,N)=M(L,N)/M(N,N)
       D(L)=M(L,L)
       U(L,L)=1.
      ENDDO

      DO L=N-1,1,-1
       DO J=L+1,N
        D(L)=D(L)-U(L,J)*U(L,J)*D(J)
       ENDDO
       IF (L.GT.1) THEN
        DO K=1,L-1
         U(K,L)=M(K,L)
         DO J=L+1,N
          U(K,L)=U(K,L)-U(K,J)*U(L,J)*D(J)
         ENDDO
         U(K,L)=U(K,L)/D(L)
        ENDDO
       ENDIF
      ENDDO

      DO K=1,N
       DO L=1,N
        V(K,L)=0.
       ENDDO
       V(K,K)=1.
      ENDDO
      
      DO K=1,N  
       DO L=N-1,1,-1
        DO J=L+1,N             
         V(L,K)=V(L,K)-U(L,J)*V(J,K)
        ENDDO
       ENDDO
      ENDDO

      DO K=1,N
      DO L=1,N
       MINV(K,L)=0.
       DO J=1,N
        MINV(K,L)=MINV(K,L)+V(J,K)*V(J,L)/D(J)
       ENDDO
      ENDDO
      ENDDO

      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 PL2U(CX,CY,CZ,CAX,IMAX,TRJ_AXIS,TRJ,
     .                POL_AXIS,POL,MMFL,K11,PX0,LN0,PX,LN,U)
c   ................................................

      IMPLICIT NONE

      INTEGER           IMAX
      
      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  MMFL
      DOUBLE PRECISION  K11
      DOUBLE PRECISION  PX0
      DOUBLE PRECISION  LN0
      DOUBLE PRECISION  PX
      DOUBLE PRECISION  LN
      DOUBLE PRECISION  U(3)
      DOUBLE PRECISION  WX(3)
      DOUBLE PRECISION  WY(3)
      DOUBLE PRECISION  WZ(3)
      DOUBLE PRECISION  Z1
      
      CALL POINTING(CX,CY,CZ,CAX,IMAX,TRJ_AXIS,TRJ,
     .              POL_AXIS,POL,LN0,LN,WX,WY,WZ)

      Z1=(PX-PX0)/K11           ! PXMM
      U(1)=Z1*WX(1)+MMFL*WZ(1)
      U(2)=Z1*WX(2)+MMFL*WZ(2)
      U(3)=Z1*WX(3)+MMFL*WZ(3)
      CALL VHAT(U,U)
          
      RETURN
      END

c  ...............................................................
      subroutine U2V(UZ,V)
c  ...............................................................

      IMPLICIT NONE
      
      real*8        w1, w2, a, b, c,
     .              x(0:3), y(0:3), z(0:3), z0, N(3), eps,
     .              p0, p1, p2, p3, z1,
     .              vec(3,0:512,0:512,6), a0(0:512,0:512,6)

      integer*4     q, i, j, f, i1, j1, l
      
      DOUBLE PRECISION      VDOT
      DOUBLE PRECISION      VNORM
      DOUBLE PRECISION      V(3)
      DOUBLE PRECISION      UX(3)      
      DOUBLE PRECISION      UY(3)
      DOUBLE PRECISION      UZ(3)
      DOUBLE PRECISION      MX2

      CHARACTER*72          INFILE

      INTEGER*4             IFF
      SAVE                  IFF, Q, VEC, MX2
      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

      eps=0.01

      Z0=1
      CALL VHAT(UZ,UZ)
      DO I=1,3
        IF(ABS(UZ(I)).LT.Z0) THEN
          J=I
          Z0=ABS(UZ(I))
        ENDIF
        UX(I)=0
      ENDDO
      UX(J)=1      
      CALL UCRSS(UZ,UX,UY)
      CALL UCRSS(UY,UZ,UX)

      Z0=-1.d10
      do f=1,6
      do i=0,q-1
      do j=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)
          if(VDOT(uz,v).lt.0) go to 10
          x(l)=VDOT(ux,v)
          y(l)=VDOT(uy,v)
          IF((x(l)**2+y(l)**2).gt.mx2) GO TO 10
          z(l)=VNORM(v)
        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((w1.ge.-eps).and.(w1.le.1+eps).and.
     .         (w2.ge.-eps).and.(w2.le.1+eps)) then
              if((z(0)+z(1)*w1+z(2)*w2+z(3)*w1*w2).gt.z0) THEN
                Z0=z(0)+z(1)*w1+z(2)*w2+z(3)*w1*w2 
              endif
            endif
          endif
        endif
10      CONTINUE
      enddo
      enddo
      enddo

      CALL VSCL(Z0,UZ,V)

      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
      character*72     infile
      character*60     line
	  
      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='(a60)') line
          if(line(43:43).eq.' ') then
            read(line,*) (vec(k,i,j,f), k=1,3)
          else
            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 FIND_PICFILE(PICNM,PICFILE,EX)
c  ...............................................................

      IMPLICIT NONE
      
      INTEGER*4         I
      INTEGER*4         SLEN
      CHARACTER*12      PICNM
      CHARACTER*72      PICFILE
      LOGICAL           EX

      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
      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
      PICFILE='NONE'

      RETURN
      END

c   ..................................................
      subroutine FLP(n,lsb,ch1,ch2)
c   ..................................................

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

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


