c
c      MARS RENDERING
c           by Dan Bruton (astro@sfasu.edu)
c
c      This program will generate an image of Mars for a given
c      viewpoint and Sun location.  Limb darkening effects are included.
c
c      This program can be found at 
c      http://www.isc.tamu.edu/~astro/dis/render/index.html
c      and was last updated on August 17, 1996.
c
c      NetPBM's ppmtogif can be used to convert the ppm image generated
c      to a gif.  The red, green and blue values (RGB) are
c      assumed to vary linearly with wavelength (for GAM=1).  
c      NetPBM Software: ftp://ftp.cs.ubc.ca/ftp/archive/netpbm/
c
       IMPLICIT REAL*8 (a-h,o-z)
       CHARACTER *8 A
       DIMENSION MP(256,128,3), ICV(700,700,3)
c
c      IMAGE INFO: WIDTH, HEIGHT, DEPTH, 
c
       SCALE=1.
       M=256*SCALE
       N=M
       MX=255
c
c      RADIUS AND POSITON OF THE EARTH IN PIXELS
c
       RAD=100.*SCALE
       XP=M/2.
       YP=N/2.
c
c      LATITUDE AND LONGITUDE OF THE OBSERVER
c      (+ east, - west)
c      AND THE TILT OF THE AXIS (DELTA)
c
       PI=3.141592654D0
       RLAT=30.*PI/180.
       RLON=-96.*PI/180.
       DELTA=15.*PI/180.
C
c      LATITUDE AND LONGITUDE OF THE SUN
c
       SLAT=30.*PI/180.
       SLON=-96.*PI/180.
c
c      READ LATITUDE AND LONGITUDE COLOR VALUES
c
       PRINT *, 'Reading data...'
       OPEN(UNIT=20,FILE='marsmap1.ppm',STATUS='UNKNOWN')
2      FORMAT (A2)
3      FORMAT (I4,I4)
4      FORMAT (I4)
5      FORMAT (I3,I4,I4)
       READ (20,2) A
       READ (20,3) II,JJ
       READ (20,4) KK
       DO J=1,128
        DO I=1,252,6
         READ (20,*) MP(I,J,1),MP(I,J,2),MP(I,J,3),
     *               MP(I+1,J,1),MP(I+1,J,2),MP(I+1,J,3),
     *               MP(I+2,J,1),MP(I+2,J,2),MP(I+2,J,3),
     *               MP(I+3,J,1),MP(I+3,J,2),MP(I+3,J,3),
     *               MP(I+4,J,1),MP(I+4,J,2),MP(I+4,J,3),
     *               MP(I+5,J,1),MP(I+5,J,2),MP(I+5,J,3)
c         PRINT *, I, MP(I,J,1),MP(I,J,2),MP(I,J,3)
        ENDDO
        READ (20,*) MP(253,J,1),MP(253,J,2),MP(253,J,3),
     *              MP(254,J,1),MP(254,J,2),MP(254,J,3),
     *              MP(255,J,1),MP(255,J,2),MP(255,J,3),
     *              MP(256,J,1),MP(256,J,2),MP(256,J,3)
       ENDDO

       DO J=1,128
        DO I=1,256
c          READ(20,5) MP(I,J,1),MP(I,J,2),MP(I,J,3)
        ENDDO
       ENDDO
       CLOSE(20)
c
c      FIND COLOR VALUE, ICV, OF EACH PIXEL
c
       PRINT *, 'Mapping...'
       RMAXVAL=0.D0
       DO J=1,N
        DO I=1,M
         DO K=1,3
          RHO=DSQRT((I-XP)**2+(J-YP)**2)
          IF(RHO.LE.RAD) THEN
            TX=I-XP
            TY=J-YP
            TZ=DSQRT(RAD**2-TX**2-TY**2)
C          ROTATION BY DELTA
            TXP=TX*DCOS(DELTA)-TY*DSIN(DELTA)
            TYP=TX*DSIN(DELTA)+TY*DCOS(DELTA)
            TZP=TZ
C          ROTATION BY RLAT 
            TXPP=TXP
            TYPP=TYP*DCOS(RLAT)-TZP*DSIN(RLAT)
            TZPP=TYP*DSIN(RLAT)+TZP*DCOS(RLAT)
C          FIND POSTION ON MAP (II,JJ)   
            THETA=DASIN(TYPP/RAD)
            JJ=1+INT((THETA/PI + 0.5)*128.)
            PHI=DASIN(TXPP/(RAD*DCOS(THETA)))+RLON
            II=1+INT((PHI/(2.*PI) + 0.5)*256.)
            IF (TZPP.LT.0.) THEN
               PHI=DASIN(-TXPP/(RAD*DCOS(THETA)))+RLON
               II=1+INT((PHI/(2.*PI) + 0.5)*256.)
               II=128+II
            ENDIF
           IF (II.GT.256) II=II-256
           IF (II.LT.1) II=II+256
            IF ((II.LE.256).AND.(II.GE.1).AND.(JJ.LE.128)
     *             .AND.(JJ.GE.1)) THEN
c
c             LIMB DARKENING
c      
              ANG=-SLON+RLON-PI
              DARKX=TXPP*DCOS(SLAT)*DSIN(ANG)/RAD
              DARKY=-TYPP*DSIN(SLAT)/RAD
              DARKZ=-TZPP*DCOS(SLAT)*DCOS(ANG)/RAD
              DARK=DARKX+DARKY+DARKZ
              IF (DARK.LT.0.D0) DARK=0.
              IF (DARK.GT.1.D0) DARK=1.
              ICV(I,J,K)=INT(DARK*MP(II,JJ,K))
            ENDIF
          ELSE
            ICV(I,J,K)=0
          ENDIF
          IF (ICV(I,J,K).GT.RMAXVAL) RMAXVAL=ICV(I,J,K)
         ENDDO
        ENDDO
       ENDDO
c
c*****************************************************************************
c
c      WRITE OUTPUT TO PPM FILE
c
c      Output color values are normalized to color depth, MX.
c
       PRINT *, 'Writing output...' 
       OPEN(UNIT=20,FILE='temp.ppm',STATUS='UNKNOWN')
1      FORMAT (A10)
       WRITE(20,1) 'P3        '
       WRITE(20,1) '# temp.ppm'
       WRITE(20,*) M,N
       WRITE(20,*) MX
       DO J=1,N
        DO I=1,M
         DO K=1,3
          ICV(I,J,K)=INT(MX*ICV(I,J,K)/RMAXVAL)
         ENDDO
         WRITE(20,*) ICV(I,J,1),ICV(I,J,2),ICV(I,J,3)
        ENDDO
       ENDDO
       CLOSE (20)
       STOP
       END 
c
c ****************************************************************************
c