c  gfortran dumpMapInc.f -O3

c  Version 1.0 - 1 Dec 2017 - Copied from dumpMapLRO_NewEmiss.f; modified to generate incident angle based on Azimuth and Zenith, does nothing else.
c Version 1.1 - 6 Dec 2017 - Added the ability to enter subsolar lat and long


      IMPLICIT NONE

      INTEGER               NTMP
      PARAMETER            (NTMP=2000)

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

      REAL*4                HT0(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                AL0(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                slope(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                inc(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                emm(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                baseEmm(-NTMP:NTMP,-NTMP:NTMP)
      REAL                  tempA (-NTMP:NTMP,-NTMP:NTMP)
      REAL                  fluxA (-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                thermSpec(30000)
      REAL*4                dHx, dHy
      REAL*4                mVal, albScale
      REAL*4                minV, delta, holdMin
      REAL*4                a, b, c, cosA, cosC, Az, Zn
      REAL*4                bigC, bigB, doubleB, bigD
      REAL*4                aPrime, bPrime, cPrime, dPrime
      REAL*8                SUNx, SUNy, SUNz, Nx, Ny, Nz
      REAL*8                temp, SUNm, Nm
      REAL*4                angle, avgTemp, avgAlb, avgInc
      REAL*4                azmuth, elevation
      REAL*4                sunLat, sunLon, objLat, objLon
      REAL*4                spaceLat, spaceLon
      REAL*4                Pi
      REAL*4                LonInput
      REAL*4                solarDist
      REAL*4                g_version
      REAL*4                targetLat, targetLon, zz, xx, yy
      REAL*4                subsolarLat, subsolarLon, t0, t1, p0, p1
      REAL*4                  phi1, phi2, tanAlpha


      INTEGER               QSZ
      INTEGER               QSZLESS
      INTEGER               I
      INTEGER               J
      INTEGER               K
      INTEGER               Istart, Iend, Jstart, Jend
      INTEGER               dX, dY, d
      INTEGER               cnt
      integer               fluxF, tempF, iF, eF, aF, slopeF
      integer               whichImage
      integer               argc
      integer               evinitX, evflipX, flip, nslope, whichway

      CHARACTER*6           MAP0
      CHARACTER*72          LMRKFILE
      CHARACTER*72          outName
      character*5          whichImgStr
      character*80          argv
      CHARACTER*72          EmmName
      g_version = 1.1
      Pi = 3.1415926

      argc = IArgC()
      if (argc .NE. 3) then
         WRITE(6,*) '# Input input map name'
         READ(5,FMT='(A6)') MAP0
         WRITE(6,*) '# Calculate based upon subsolar (1) or az-zn (!1)?'
         READ(5,*) whichway
         if (whichway .eq. 1) then
          WRITE(6,*) '# Input subsolar lat'
          READ(5,*) sunLat
          WRITE(6,*) '# Input subsolar lon'
          READ(5,*) sunLon
          WRITE(6,*) '# Input lat of center of map'
          READ(5,*) objLat
          WRITE(6,*) '# Input lon of center of map'
          READ(5,*) objLon
         else 
          WRITE(6,*) '# Azimuth (in deg)'
          READ(5,*) Az
          WRITE(6,*) '# Zenith (in deg)'
          READ(5,*) Zn
         endif
      else
         call getArg (1, MAP0)
         call getArg (2, argv)
         READ(argv,*) whichway
         if (whichway .eq. 1) then
          call getArg (3, argv)
          READ(argv,*) sunLat
          call getArg (4, argv)
          READ(argv,*) sunLon
          call getArg (5, argv)
          READ(argv,*) objLat
          call getArg (6, argv)
          READ(argv,*) objLon
         else
          call getArg (3, argv)
          READ(argv, *) Az
          call getArg (4, argv)
          READ(argv, *) Zn
         endif
      endif

      if (whichway .eq. 1) then
       t0 = (90 - objLat) * Pi/180
       t1 = (90 - sunLat) * Pi/180
       p0 = objLon * Pi/180
       p1 = sunLon * Pi/180
       zz = cos(t0)*cos(t1) + sin(t0)*sin(t1)*cos(p1 - p0) 
       xx = sin(t1)*sin(p1 - p0)
       yy = sin(t0)*cos(t1) - cos(t0)*sin(t1)*cos(p1 - p0)

       Az = ATAN(xx/yy) * 180/Pi
       Zn = ACOS(zz) * 180/Pi
      endif

      WRITE(6,*) 'MAP0, Az, Zn: ', MAP0, Az, Zn

      LMRKFILE='MAPFILES/'//MAP0//'.MAP'
      CALL READ_MAP(LMRKFILE,NTMP,QSZ,SCALE,V,UX,UY,UZ,HT0,AL0)

      WRITE(6,*) '# QSize ', QSZ
      if (qsz .gt. NTMP) then
        write (*,*) "Qsize is too big for this program: ", qsz, NTMP
        stop
      endif

c     Convert Az and Zn to radians
      Az = Az * Pi / 180
      Zn = Zn * Pi / 180

      WRITE(6,*) 'Az Zn: ', Az, Zn


c ########################################
c Slope
c ########################################
c		 Loop to sample a large area around what I want to see

c     Wide field of view
      Istart = -qsz
      Iend = qsz
      Jstart = -qsz
      Jend = qsz

      WRITE(6,*) '# Version: ', g_version
      WRITE(6,*) '# Calculating slope'
      DO J=Jstart,Jend
      DO I=Istart,Iend

c     Incident angle
c     Calculate sun vector

      SUNx = sin (Az) * sin (Zn)
      SUNy = cos (Az) * sin (Zn)
      SUNz = cos (Zn)
      SUNm = SQRT(SUNx*SUNx + SUNy*SUNy + SUNz*SUNz)

c     Calculate normal at pixel
c     Bring back dHx and dHy

c     Change HT0 values so we can see the math easily
c
c      HT0(I+1,J)=2
c      HT0(I-1,J)=2
c      HT0(I,J+1)=2
c      HT0(I,J-1)=0
c
c     END Changing of HT0 values

      dHx = (HT0 (I+1,J) - HT0 (I-1,J) )
      dHy = -(HT0 (I, J+1) - HT0 (I, J-1) )


c     Cross product results from y slope vector and x slope vector
      Nx = -dHx * 2.0
      Ny = -dHy * 2.0
      Nz = 4.0
      Nm = SQRT(Nx*Nx + Ny*Ny + Nz*Nz)

c     Testing only below take out before final analysis!
      if (I .eq. 199) then
       if (J .eq. 461) then
        inc(I,J) = ACOS((Nx * SUNx + Ny * SUNy + Nz * SUNz) / (Nm*SUNm))
        WRITE(6,*) 'Normal Vector', Nx, Ny, Nz, Nm
        WRITE(6,*) 'Normal Vector Sq', Nx*Nx, Ny*Ny, Nz*Nz
        WRITE(6,*) 'HTO J direction', HT0(I, J+1), HT0(I, J-1)
        WRITE(6,*) 'HTO I direction', HT0(I+1,J), HT0(I-1,J)

        WRITE(6,*) 'Sun vector (x,y,z) ', SUNx, SUNy, SUNz, SUNm
        WRITE(6,*) 'Incident Angle vec x,y,z', Nx*SUNx, Ny*SUNy, Nz*SUNz
        WRITE(6,*) 'Incident Angle', inc(I,J)*180/Pi

      endif
      endif

c      Nx = 0
c      Ny = Nz  
c      Nm = SQRT(Nx*Nx + Ny*Ny + Nz*Nz)
c     End Testing

c     normal vector dotted with Sun vector
      inc(I,J) = ACOS((Nx * SUNx + Ny * SUNy + Nz * SUNz) / (Nm*SUNm))

c     Solar vector and object vector
c      write (*,*) "Sun (C): ", sunLat*180/Pi, sunLon*180/Pi, 
c     +      "b", C*180/Pi
c      write (*,*) "Object (B): ", objLat*180/Pi, objLon*180/Pi, 
c     +      "c", B*180/Pi


c      write (*,*) "cosA", cosA
c      write (*,*) "cosA angle ", acos(cosA)*180/Pi

      ENDDO
      ENDDO

c ########################################
c Print Grid
c ########################################
      DO J=Jstart,Jend
        DO I=Istart,Iend
          inc(I,J) = inc(I,J) * 180 / 3.1415926

        ENDDO
      ENDDO
      outName=MAP0//'-i.txt'
      CALL write_array (inc, NTMP,QSZ, outName)

      STOP
      END

c   ................................................
      SUBROUTINE planck(currTemp, spec)
c   ................................................
      IMPLICIT NONE
      
      REAL               currTemp
      REAL*4             spec(30000)
      INTEGER            wave
      REAL*4             h, c, boltz, flux, emm, st_boltz, tWave
      REAL*4             top, bottom, eTerm

      h = 6.626069e-34
      c = 299792458
      boltz = 1.380648e-23
      flux = 1370
      emm = 0.94
      st_boltz = 5.67e-8

      DO wave=1,30000
        tWave = wave/1.0e9
        top = 2.0 * h * c * c
        eTerm = h * c / tWave / boltz / currTemp
        bottom = tWave**5 * (exp (eTerm) - 1)
       
        flux = top/bottom
        spec(wave) = spec(wave) + flux

c        write (*,*) wave/1000.0, flux
      enddo

  
      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')

c       Read is the scale fo the bigmap
        READ(20,REC=1) BLINE
        CH4f=BLINE(7:10)
        call flip(4,lflag,ch4f,ch4)
        SCALE=RL4

c       Read is the size of the bigmap
        QSZ=ICHAR(BLINE(11:11))
     .     +ICHAR(BLINE(12:12))*256

c       Loop over K four times (for values for V, and U[XYZ])
c       Load up V, UX, UY and UZ
        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


c       Load the H Scale
        CH4f=BLINE(64:67)
        call flip(4,lflag,ch4f,ch4)
        HSCALE=RL4

c       Clear the memory
        DO I=-QSZ,QSZ
          DO J=-QSZ,QSZ
            ALB(I,J)=0
            HT(I,J)=0
          ENDDO
        ENDDO
        
c       Read stuff in, 24 bytes at a time
        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

c       Read the remainder of a 24 byte line
        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


c   ................................................
      SUBROUTINE write_array(temp0, NTMP, QSZ, fName)

c   ................................................

      IMPLICIT NONE
      
      INTEGER               NTMP

      INTEGER               NREC
      INTEGER               LENSIZE
      DOUBLE PRECISION      SCALE
      DOUBLE PRECISION      HSCALE
      DOUBLE PRECISION      V(3)
      INTEGER               QSZ
      INTEGER               I
      INTEGER               J
      REAL                temp0(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                sunLon


      CHARACTER*72          fName

      WRITE(6,*) '# Output: ', fName
      OPEN(UNIT=10, FILE=fName)
      DO J=-qsz,qsz
        DO I=-qsz,qsz
          write (10,"(f15.4)", advance='no') temp0(I,J)
        enddo
        write (10,*)
      enddo
      CLOSE(UNIT=10)

      RETURN
      END     

c   ................................................
      SUBROUTINE read_array(temp0, NTMP, QSZ, fName)

c   ................................................

      IMPLICIT NONE

      INTEGER               NTMP

      INTEGER               NREC
      INTEGER               LENSIZE
      DOUBLE PRECISION      SCALE
      DOUBLE PRECISION      HSCALE
      DOUBLE PRECISION      V(3)
      INTEGER               QSZ
      INTEGER               I
      INTEGER               J
      REAL                temp0(-NTMP:NTMP,-NTMP:NTMP)
      REAL*4                sunLon


      CHARACTER*72          fName

      OPEN(UNIT=10, FILE=fName)
      DO J=-qsz,qsz
        DO I=-qsz,qsz
          write (10,"(f15.4)", advance='no') temp0(I,J)
        enddo
        write (10,*)
      enddo
      CLOSE(UNIT=10)

      RETURN
      END


 
