; MAIN PROGRAM
; contour plots of duct flow for the program duct.
; this assumes the grid is cartesian and uniform
; the program is derived from fivol plotting - thus some funny names

   name='' & contin='' & again='y' & withps='n' & run='' & fnumber=1
   pi=3.141593

; GRID 
   nx = long(51)  &  ny = long(51) &  niter = long(51)  &  me=long(1) 
   flowc=0.1 &  flowx=0.1  &  rms=0.1 

; READ PARAMETERS 
   openr, 8, 'duct.bin',/F77_UNFORMATTED
   readu, 8,  nx,ny
   print, 'dimension nx=',nx,'     ny=',ny
   phi=fltarr(nx,ny) & x=fltarr(nx) & y=fltarr(ny) & phix=phi & source=phi 
   readu, 8,  x,y
   print, ' x, y read'
   readu, 8,  phix,phi
   print, ' phix,phi read'
   readu, 8,  flowc,flowx
   print, ' flowc,flowx read'
   readu, 8,  me,niter,rms,omega
   print, ' niter, rms om read'
   close,8
   

    print, 'Which case?'
    read, run

     while (again eq 'y') do begin

      !P.THICK=1.
; POSTSCIPT OUTPUT      
      print, 'With postscript?'
      read, withps
      if withps eq 'y' then begin 
        set_plot,'ps'
        device,filename='con.ps'
        !P.THICK=2.
        device,/inches,xsize=8.,scale_factor=1.0,xoffset=0.5
        device,/inches,ysize=10.0,scale_factor=1.0,yoffset=0.5
        device,/times,/bold,font_index=3
       endif

; COORDINATES OF PLOT AND WRITTEN OUTPUT
        xmax=max(x) & xmin=min(x) & dxpos=(xmax-xmin)
        ymax=max(y) & ymin=min(y) & dypos=(ymax-ymin)
       
        xpos=xmax + 0.05*dxpos
        ypos=ymin + 0.96*dypos  
        ypos0=ymin + 0.9*dypos  
        ypos1=ymin + 0.84*dypos   
        ypos2=ymin + 0.78*dypos   
        ypos3=ymin + 0.72*dypos   
        ypos4=ymin + 0.6*dypos   
        ypos5=ymin + 0.54*dypos   
        ypos6=ymin + 0.48*dypos   
        ypos7=ymin + 0.42*dypos   
         
        xa1 = 0.1 & xe1 = 0.6
        ylo1 = 0.5 & yup1 = 0.9
        ylo2 = 0.07 & yup2 = 0.47

       
  print, 'plot  page? '
  read, contin
  if (contin eq '' or contin eq 'y') then begin

; PLOT PAGE
       !P.REGION=[0.,0.,1.0,1.0]
       !P.MULTI=[0,2,0,0,0]
       !P.CHARSIZE=1.0
       !P.FONT=3
;       !X.TICKS=2
;       !Y.TICKS=8
;       !Y.TICKlen=0.04
       !X.THICK=2
       !Y.THICK=2
       !P.THICK=2.
       !X.RANGE=[xmin,xmax]
       !Y.RANGE=[ymin,ymax]

        nlevels=11  & lvec=findgen(nlevels) 
        dlvec=(max(phix)-min(phix))/(nlevels-1)
        lvec=lvec*dlvec+min(phix) 
        maxphi=max(phi) & maxphix=max(phix)
        
	!P.POSITION=[xa1,ylo1,xe1,yup1]
	contour,phix,x,y,levels=lvec,$
        title=' Flow W (Exact solution dotted)',$
        xstyle=1,ystyle=1,c_linestyle=1
	contour,phi,x,y,levels=lvec,xstyle=1,ystyle=1
        if (me eq 1) then xyouts,xpos,ypos,'Method: FEM' $
             else xyouts,xpos,ypos,'Method: FDM' 
	xyouts,xpos,ypos0,$
	    'nx='+string(nx,'(i2)')+'   ny='+string(ny,'(i2)')
	xyouts,xpos,ypos1,'Omega: '+string(omega,'(f6.2)')
	xyouts,xpos,ypos2,'Error: '+string(rms,'(f7.5)')
	xyouts,xpos,ypos3,'No. iter.: '+string(niter,'(i3)')
	xyouts,xpos,ypos4,'Max computed flow: '+string(maxphi,'(f7.5)')
	xyouts,xpos,ypos5,'Max analytic flow: '+string(maxphix,'(f7.5)')
	xyouts,xpos,ypos6,'Computed flowrate: '+string(flowc,'(f7.4)')
	xyouts,xpos,ypos7,'Analytic flowrate: '+string(flowx,'(f7.4)')

;        PLOT, x, y, psym=1,xstyle=4,ystyle=4
;        axis, xax=0, -0.05*dxpos,-0.1*dypos, xtitle='X'
;        axis, yax=0, -0.05*dxpos,-0.1*dypos, ytitle='Y'

        !P.THICK=3.
;        plots, bx1,by1,/data
;        plots, bx2,by2,/data
;        plots, bx3,by3,/data
;        plots, bx4,by4,/data
        
        
  endif

     !P.FONT=3


     print, 'view results again or make ps file?'
     read, again
     if withps eq 'y' then device,/close
     set_plot,'x'
     !P.THICK=1.
   endwhile

end