PHOTON USE
  p;parphi
  1 2.5 1 ;;;;
 
  msg The grid. The z coordinate represents time
  msg Press return for the pressure contours.
  gr x 1;pause
  gr off;red;gr ou x 1
  msg pressure distribution. Press return for enthalpies
  con p1 x 1 fi;0.1;pause;con off;red
  msg enthalpy distribution. Press e to end, or type other command
  con h1 x 1 fi;0.1
  enduse
TEXT(1D Gas Oscillation; Various Cases  
TITLE
  DISPLAY
 
  The default settings are for a pipe with closed ends, and with
  initially uniform-pressure and -temperature gas in uniform motion.
 
  A series of intersecting pressure waves results, with
  corresponding variations of enthalpy.
 
  Other cases involve the gas initially at rest, and the presence
  of time-varying boundary conditions at one end of the pipe. In one
  case the pressure is fixed at an intemediate point.
 
  Ideal-gas properties are used.
 
  PHOTON USE commands are provided.
 
  ENDDIS
char(casename)
real(vin,win,tin,pin,gasconst,hin,cp,cv,gamma)
gasconst=1.0;gamma=1.4;cv=gasconst/(gamma-1.0);cp=cv*gamma
vin=0.5;win=1.0;tin=1.0;hin=cp*tin;pin=0.0;lstep=500
delay(100)
mesgm(This case OK? If not, press n
readvdu(ans,char,y)
integer(kase)
kase=0
if(:ans:.eq.n) then
 mesga(Please choose from the following:-
 mesg (1. Gas at rest. Sudden sustained entry from top.
 mesg (2. Gas at rest. Intermittent entry from top.
 mesg (3. As for 2, with pressure=press0 at ny/4
 readvdu(kase,int,1)
 mesg(Case number :kase: has been chosen
endif
 
if(kase.gt.0) then
 vin=0.0
endif
steady=f
grdpwr(t,lstep,5.0,1.0)
 
   GROUP 4. Y-direction grid specification
GRDPWR(Y,20,1.0,1.0)
LSTEP
NY
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,h1)
 
    GROUP 9. Properties of the medium (or media)
RHO1=IDEALGAS;RHO1A=0.0;RHO1B=1.0;PRESS0=1.0;RHO1C=0.0
DRH1DP=RHO1B;TMP1=LINH;TMP1B=1.0/CP; CP1=CP
 
    GROUP 11. Initialization of variable or porosity fields
FIINIT(P1)=PIN;FIINIT(V1)=VIN;FIINIT(H1)=HIN
 
    GROUP 13. Boundary conditions and special sources
 
if(kase.gt.0) then
 real(floin)
 floin=0.5
 patch(timn,north,1,1,ny,ny,1,1,1,lstep)
 coval(timn,p1,fixflu,grnd1)
 coval(timn,v1,onlyms,-floin)
 coval(timn,h1,onlyms,4*hin)
 itima=2*lstep
 tima=floin
 casename=Gas_at_rest._Sudden_sustained_entry_from_top
endif
 
if(kase.gt.1) then
 floin=0.1
 patch(timn,north,1,1,ny,ny,1,1,1,lstep)
 mesg(Intermittent gas entry and withdrawal at top of pipe
 mesg(according to...
 mesgb(1. battlement ;  2. saw-tooth;  3. sine curve ? Insert number
 integer(number);real(arg)
 arg=grnd1
 readvdu(number,int,1)
 if(number.eq.2) then
+ arg=grnd2
 endif
 if(number.eq.3) then
+ arg=grnd3
 endif
 coval(timn,p1,fixflu,arg)
 coval(timn,v1,onlyms,-floin)
 coval(timn,h1,onlyms,4*hin)
 tima=floin
 casename=Gas_at_rest._Intermittent_entry_from_top
 itima=lstep/4;itimc=lstep
endif
 
if(kase.gt.2) then
 casename=pressure_fixed_at_quarter_height
 patch(pfix,cell,1,1,ny/4,ny/4,1,nz,1,lstep)
 coval(pfix,p1,fixval,0.0)
 coval(pfix,v1,onlyms,0.0)
 coval(pfix,h1,onlyms,hin)
endif
 
    GROUP 16. Termination of iterations
LSWEEP=10
SELREF=T;RESFAC=1.E-3
VARMIN(H1)=0.5
 
    GROUP 21. Print-out of variables
 
    GROUP 22. Spot-value print-out
IXMON=1;IYMON=1
 
    GROUP 23. Field print-out and plot control
TSTSWP=LSWEEP
ITABL=1;NYPRIN=1;NTPRIN=LSTEP/4
IDISPA=10;IDISPB=1;IDISPC=LSTEP
PATCH(PROFILES,PROFIL,1,1,NY/2,NY/2,1,1,1,LSTEP)
PLOT(PROFILES,P1,0.0,0.0)
PLOT(PROFILES,V1,0.0,0.0)
text(:casename: