C     PROGRAM: The Phong reflection model.
C
C     Dolocitev spremenljivk.
      INTEGER nX,nY,Width,Height,Square,Radius,ChoiceLight
      INTEGER x,y,RGB,Pic(1024,1024,3),R,G,B,Rback,Gback,Bback
      INTEGER Xdelta,Ydelta,Xtable(10,10),Ytable(10,10),Xc,Yc
      INTEGER SpecIndex,WrM(3145728),ColourMax,ChoiceColour,ChoiceBack
      REAL    Kd,Ks,Kdelta,Step,Pi,Fi,Theta,Lx,Ly,Lz
C     Konstante
      PARAMETER(Pi=3.141593)
C     Odpiranje datoteke Phong.ppm  
      OPEN(1,FILE='Phong.ppm')
C     Naslovna stran programa PHONG.
      WRITE(*,*)
      WRITE(*,*)
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'*               ~ Program  PHONG.EXE ~               *'
      WRITE(*,*)'*          >> The Phong reflection model <<          *'
      WRITE(*,*)'*                Avtor: Bostjan VEBER                *'
      WRITE(*,*)'*    --------------------------------------------    *'
      WRITE(*,*)'*    Serija osencenih krogel s svetlobnim modelom    *'
      WRITE(*,*)'*     Phong pri razlicnih vrednostih parametrov:     *'
      WRITE(*,*)'*                                                    *'
      WRITE(*,*)'*       # Kd - parameter razprsitve svetlobe         *'
      WRITE(*,*)'*       # n  - parameter sirine odbleska             *'
      WRITE(*,*)'*       # Ks - parameter moci odbleska               *'
      WRITE(*,*)'******************************************************'
      WRITE(*,*)
C
C     Dolocitev stevila krogel v X smeri.
10    WRITE(*,*)
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'[1] Podaj stevilo krogel v X smeri (max 10 krogel):   '
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)'nx =?'
      READ(*,*)nX
      IF((nX.GT.10).OR.(nX.LT.1))THEN
      GOTO 10
      END IF
C
C     Dolocitev stevila krogel v Y smeri.
20    WRITE(*,*)
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'[2] Podaj stevilo krogel v Y smeri (max 10 krogel):   '
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)'ny =?'
      READ(*,*)nY
      IF((nY.GT.10).OR.(nY.LT.1))THEN
      GOTO 20
      END IF
C
C     Dolocitev sirine risbe (nX>=nY).
      IF(nX.GE.nY)THEN
30    WRITE(*,*)
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'[3] Podaj sirino risbe v pt (min 320, max 1024):      '
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)'Sirina =?'                             
      READ(*,*)Width
      IF((Width.GT.1024).OR.(Width.LT.320))THEN
      GOTO 30
      END IF      
C     Izracun osnovnih parametrov risbe.
      Square=Width/nX
      Height=Square*nY
      Radius=Square*0.4
C      
C     Dolocitev visine risbe (nY>nX).
      ELSE
40    WRITE(*,*)
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'[3] Podaj visino risbe v pt (min 320, max 1024):      '
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)'Visina =?'
      READ(*,*)Height
      IF((Height.GT.1024).OR.(Height.LT.320))THEN
      GOTO 40
      END IF
C     Izracun osnovnih parametrov risbe.
      Square=Height/nY
      Width=Square*nX
      Radius=Square*0.4
      END IF
C      
C    Dolocitev parametra razprsitve svetlobe Kd pri (nX=1).
      IF(nX.EQ.1)THEN
50    WRITE(*,*)
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'[4] Podaj vrednost parametra razprsitve svetlobe <Kd> '
      WRITE(*,*)'    med 0 in 1. Vrednost parametra parametra moci     '
      WRITE(*,*)'    odbleska <Ks> je dolocena z enacbo Ks=1-Kd.       '
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)'Kd =?'
      READ(*,*)Kd
      IF((Kd.GT.1).OR.(Kd.LT.0))THEN
      GOTO 50
      END IF
      Ks=1-Kd
      Kdelta=0
C     Dolocitev parametra razprsitve svetlobe Kd pri (nX>1).
      ELSE
      Kd=0
      Ks=1
      Kdelta=1.0/nX
      END IF
C
C     Dolocitev parametra sirine odbleska n pri (nY=1).     
      IF(nY.EQ.1)THEN
60    WRITE(*,*)
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'[5] Podaj celostevilsko vrednost parametra sirine     '
      WRITE(*,*)'    odbleska <n> med 5 in 40.                         '  
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)'n =?'
      READ(*,*)SpecIndex
      IF((SpecIndex.GT.40).OR.(SpecIndex.LT.5))THEN
      GOTO 60
      Step=1
      END IF
C     Dolocitev parametra sirine odbleska n pri (nY>1).     
      ELSE
      SpecIndex=5
      Step=8.0**(1.0/(nY-1))
      END IF
C      
C     Dolocitev barve krogel.
70    WRITE(*,*)
      WRITE(*,*)'******************************************************'          
      WRITE(*,*)'[6] Izberi barvo krogel:                              '
      WRITE(*,*)'    <1> Rdece krogle                                  '
      WRITE(*,*)'    <2> Zelene krogle                                 '
      WRITE(*,*)'    <3> Modre krogle                                  '
      WRITE(*,*)'    <4> Barva krogel po lastni izbiri (RGB model)     '
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)'Izbira barve (1/2/3/4) =?'
      READ(*,*)ChoiceColour
C      
      IF(ChoiceColour.EQ.1)THEN
      R=255
      G=0
      B=0
      ELSE
      IF(ChoiceColour.EQ.2)THEN
      R=0
      G=255
      B=0
      ELSE
      IF(ChoiceColour.EQ.3)THEN
      R=0
      G=0
      B=255
      ELSE
      IF(ChoiceColour.EQ.4)THEN
C     Izbira barve krogel po lastni izbiri.
      WRITE(*,*)
      WRITE(*,*)'******************************************************'          
      WRITE(*,*)'             >> Dolocitev barve KROGEL <<             '
      WRITE(*,*)'[7] Podaj barvo krogel v RGB modelu (Red/Green/Blue). '
      WRITE(*,*)'    Stevilke osnovnih treh barv so med 0 in 255.      '
      WRITE(*,*)'------------------------------------------------------'
C     Rdeca barva - Red.
80    WRITE(*,*)
      WRITE(*,*)'RED =?'
      READ(*,*)R
      IF((R.GT.255).OR.(R.LT.0))THEN
      GOTO 80
      END IF
C     Zelena barva - Green.
90    WRITE(*,*)
      WRITE(*,*)'GREEN =?'
      READ(*,*)G
      IF((G.GT.255).OR.(G.LT.0))THEN
      GOTO 90
      END IF
C     Modra Barva - Blue.
100   WRITE(*,*)
      WRITE(*,*)'BLUE =?'
      READ(*,*)B
      IF((B.GT.255).OR.(B.LT.0))THEN
      GOTO 100
      END IF
      ELSE
      GOTO 70
      END IF
      END IF
      END IF
      END IF
C
C     Dolocitev ozadja slike osencenih krogel.
110   WRITE(*,*)
      WRITE(*,*)'******************************************************'          
      WRITE(*,*)'[8] Izberi barvo ozadja risbe:                        '
      WRITE(*,*)'    <1> Belo ozadje                                   '
      WRITE(*,*)'    <2> Crno ozadje                                   '
      WRITE(*,*)'    <3> Barva ozadja po lastni izbiri (RGB model)     '
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)'Izbira ozadja (1/2/3) =?'
      READ(*,*)ChoiceBack
C      
      IF(ChoiceBack.EQ.1)THEN
      Rback=255
      Gback=255
      Bback=255
      ELSE
      IF(ChoiceBack.EQ.2)THEN
      Rback=0
      Gback=0
      Bback=0
      ELSE
      IF(ChoiceBack.EQ.3)THEN
      WRITE(*,*)
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'             >> Dolocitev barve OZADJA <<             '
      WRITE(*,*)'[9] Podaj barvo OZADJA v RGB modelu (red/green/blue). '
      WRITE(*,*)'    Stevilke osnovnih treh barv so med 0 in 255.      '
      WRITE(*,*)'------------------------------------------------------'
C     Rdeca barva - Red.
120   WRITE(*,*)
      WRITE(*,*)'RED =?'
      READ(*,*)Rback
      IF((Rback.GT.255).OR.(Rback.LT.0))THEN
      GOTO 120
      END IF
C     Zelena barva - Green.
130   WRITE(*,*)
      WRITE(*,*)'GREEN =?'
      READ(*,*)Gback
      IF((Gback.GT.255).OR.(Gback.LT.0))THEN
      GOTO 130
      END IF
C     Modra barva - Blue.
140   WRITE(*,*)
      WRITE(*,*)'BLUE =?'
      READ(*,*)Bback
      IF((Bback.GT.255).OR.(Bback.LT.0))THEN
      GOTO 140
      END IF
      ELSE
      GOTO 110
      END IF
      END IF
      END IF
C            
      DO RGB=1,3      
C     Red Colour (R).
      IF(RGB.EQ.1)THEN
      DO y=1,Height
      DO x=1,Width
      Pic(x,y,RGB)=Rback
      END DO
      END DO
      END IF
C     Green Colour (G).
      IF(RGB.EQ.2)THEN
      DO y=1,Height
      DO x=1,Width
      Pic(x,y,RGB)=Gback
      END DO
      END DO
      END IF
C     Blue Colour (B).
      IF(RGB.EQ.3)THEN
      DO y=1,Height
      DO x=1,Width
      Pic(x,y,RGB)=Bback
      END DO
      END DO
      END IF
      END DO
C     
C     Dolocitev vektorja luci.
150   WRITE(*,*)
      WRITE(*,*)'******************************************************'          
      WRITE(*,*)'[10] Izberi smer osvetlitve krogel:                   '
      WRITE(*,*)'     <1> Standardna nastavitev                        '
      WRITE(*,*)'     <2> Smer osvetlitve po lastni izbiri             '
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)'Izbira smeri osvetlitve (1/2) =?'
      READ(*,*)ChoiceLight
C      
      IF(ChoiceLight.EQ.1)THEN
      Lx=0.57735
      Ly=-0.57735
      Lz=0.57735
      ELSE
      IF(ChoiceLight.EQ.2)THEN
      WRITE(*,*)
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'           >> Dolocitev SMERI OSVETLITVE <<           '
      WRITE(*,*)'[11] Smer osvetlitve se podaja z vektorjem luci,ki ga '
      WRITE(*,*)'     dolocata kota <Fi> in <Theta> v stopinjah.       '
      WRITE(*,*)'  --------------------------------------------------  '
      WRITE(*,*)'     Lego luci doloca kot <Fi> med 0 in 360:          '
      WRITE(*,*)'       # Vzhod -> Fi=0                                '
      WRITE(*,*)'       # Sever -> Fi=90                               '
      WRITE(*,*)'       # Zahod -> Fi=180                              '
      WRITE(*,*)'       # Jug   -> Fi=270                              '
      WRITE(*,*)'                                                      '
      WRITE(*,*)'     Visino luci doloca kot <Theta> med 0 in 180:     '
      WRITE(*,*)'       # Luc je v visini ravnine XY -> Theta=0        '
      WRITE(*,*)'       # Luc je nad kroglami        -> Theta=90       '
      WRITE(*,*)'------------------------------------------------------'
      WRITE(*,*)
      WRITE(*,*)'Fi =?'
      READ(*,*)Fi
160   WRITE(*,*)
      WRITE(*,*)'Theta =?'
      READ(*,*)Theta
      IF(Theta.GT.180)THEN
      GOTO 160
      End if
      Fi=Fi*Pi/180
      Theta=Pi/2-Theta*Pi/180
      Lx=SIN(Theta)*COS(Fi)
      Ly=-SIN(Theta)*SIN(Fi)
      Lz=COS(Theta)
      ELSE
      GOTO 150
      END IF
      END IF
C
C     GLAVNI PROGRAM.
      Ydelta=Square/2
C     Rows.      
      DO I=1,nY
      Xdelta=Square/2
C     Columns.      
      DO J=1,nX
      Xtable(I,J)=Xdelta
      Ytable(I,J)=Ydelta
      Xc=Xtable(I,J)
      Yc=Ytable(I,J)
      CALL Shade(Kd,Ks,SpecIndex,Xc,Yc,Radius,Lx,Ly,Lz,Pic,R,G,B)
      Kd=Kd+Kdelta
      Ks=Ks-Kdelta
      Xdelta=Xdelta+Square
      END DO
      Ydelta=Ydelta+Square
      SpecIndex=SpecIndex*Step
      Kd=0
      Ks=1
      END DO
C      
C     Prepis matrike Pic v WrM.
      I=1
      DO y=1,Height
      DO x=1,Width
      DO RGB=1,3
      WrM(I)=Pic(x,y,RGB)
      I=I+1
      END DO
      END DO
      END DO
C      
C     Dolocitev maksimalne vrednosti indeksa barve.
      N=Width*Height*3
      ColourMax=0    
      DO I=1,N
      IF(ColourMax.EQ.255)THEN
      GOTO 200
      ELSE
      IF(WrM(I).GT.ColourMax)THEN
      ColourMax=WrM(I)
      END IF
      END IF
      END DO
C      
C     Zapis PPM - Portable Pixmap Format datoteke.
200   WRITE(1,210)
210   FORMAT('P3')
      WRITE(1,220)
220   FORMAT('# Phong.ppm')
      WRITE(1,230)Width,Height
230   FORMAT(2I5)
      WRITE(1,240)ColourMax
240   FORMAT(1I4)
      DO I=1,N-5,6
      WRITE(1,250)WrM(I),WrM(I+1),WrM(I+2),WrM(I+3),WrM(I+4),WrM(I+5)
250   FORMAT(6I4)
      END DO
C
C     Konec programa Phong.
      WRITE(*,*)
      WRITE(*,*)'******************************************************' 
      WRITE(*,*)'*    >> Program PHONG.EXE je zakljucil z delom! <<   *'
      WRITE(*,*)'*  ------------------------------------------------- *' 
      WRITE(*,*)'*  Datoteka PHONG.PPM je v istem direktoriju, kot je *'
      WRITE(*,*)'*  program PHONG.Datoteko odprete v okolju X-Windows *'           
      WRITE(*,*)'*  s programom PAINT SHOP PRO ali pa v okolju Unix s *'
      WRITE(*,*)'*  programom XV.                                     *'
      WRITE(*,*)'*  ------------------------------------------------- *'
      WRITE(*,*)'*               >> HVALA ZA UPORABO <<               *'
      WRITE(*,*)'******************************************************'
      END
C      
C     PODprogram za osencenje krogle SHADE.
      SUBROUTINE Shade(Kd,Ks,SpecIndex,Xc,Yc,Radius,Lx,Ly,Lz,Pic,R,G,B)
C
C     Spremenljivke.
      INTEGER SpecIndex,Xc,Yc,Radius,Ir,Ig,Ib,x,y,z,R,G,B
      INTEGER Pic(1024,1024,3)
      REAL Kd,Ks,rsquare,xsquare,ysquare,zsquare,denom,xn,yn,zn
      REAL LdotN,NH,NnH,Lx,Ly,Lz,dx,dy,dz,HvX,HvY,HvZ,deH,Hx,Hy,Hz
      REAL dist,distfactor,difuseterm,specularterm,D
C
C     Konstante.
      PARAMETER(Ilight=150)
      PARAMETER(K=70.0)
      PARAMETER(ambientterm=0.3)
C
C     Dolocitev vektorja oddaljenosti.
      D=Radius+130
      dx=Lx*D
      dy=Ly*D
      dz=Lz*D
C     Dolocitev vektorja H.      
      HvX=0.5*Lx
      HvY=0.5*Ly
      HvZ=0.5*(1+Lz)
      deH=SQRT(HvX**2+HvY**2+HvZ**2)
      Hx=HvX/deH
      Hy=HvY/deH
      Hz=HvZ/deH
C      
C     Glavni del PODprograma.
      rsquare=Radius**2
      DO y=-Radius,Radius,1
      ysquare=y**2
      DO x=-Radius,Radius,1
      xsquare=x**2
C      
      IF((xsquare+ysquare).LE.rsquare)THEN
      z=SQRT(rsquare-xsquare-ysquare)
      zsquare=z**2
      denom=SQRT(xsquare+ysquare+zsquare)
      xn=x/denom
      yn=y/denom
      zn=z/denom
C      
      LdotN=xn*Lx+yn*Ly+zn*Lz
      IF(LdotN.LE.0)THEN
      LdotN=0
      ELSE
      dist=SQRT((dx-x)**2+(dy-y)**2+(dz-z)**2)
      NH=Hx*xn+Hy*yn+Hz*zn
      NnH=EXP(SpecIndex*Log(NH))
      END IF
C
C     Tocke niso osvetljene od luci.
      IF(LdotN.LE.0)THEN
      Ir=R*ambientterm
      Ig=G*ambientterm
      Ib=B*ambientterm
C
C     Osvetljene tocke na krogli.
      ELSE
C      
      distfactor=Ilight/(dist+K)
      difuseterm=distfactor*Kd*LdotN
      specularterm=distfactor*Ks*NnH
C
C     Izracun osvetljenosti rdece barve - RED.      
      IF(R.EQ.0)THEN
      Ir=255*specularterm
      ELSE
      Ir=R*(ambientterm+difuseterm+specularterm)
      IF(Ir.GT.255)THEN
      Ir=255
      END IF
      END IF
C     Izracun osvetljenosti zelene barve - GREEN.      
      IF(G.EQ.0)THEN
      Ig=255*specularterm
      ELSE
      Ig=G*(ambientterm+difuseterm+specularterm)
      IF(Ig.GT.255)THEN
      Ig=255
      END IF
      END IF
C     Izracun osvetljenosti modre barve - BLUE.      
      IF(B.EQ.0)THEN
      Ib=255*specularterm
      ELSE
      Ib=B*(ambientterm+difuseterm+specularterm)
      IF(Ib.GT.255)THEN
      Ib=255
      END IF
      END IF
      END IF
C
C     Zapis matrike risbe Pic.      
      Pic(Xc+x,Yc+y,1)=Ir
      Pic(Xc+x,Yc+y,2)=Ig
      Pic(Xc+x,Yc+y,3)=Ib  
      END IF
      END DO
      END DO
C    
      RETURN
      END