inter1.f
c<html>
c<body>
c<pre>
program interpolate
implicit none
c
c Demonstation of 2 simple interpolation methods for
c smoothly joining results from 2 disjoint regions
c
c John Mahaffy 2/3/96
c
c Open up two files in your current directory to receive results
c
open(10,file='linint')
open(11,file='cubeint')
c
c First a linear interpolation
c
call linenrgy
c
c Next a cubic interpolation
c
call cubenrgy
c
c If you are in Hammond Lab, change this subroutine to plot results
c
call plotit
c
stop
end
subroutine linenrgy
c
c Generate internal energy and Cv of a gas for Temperatures
c between 300 and 3000K in steps of 10K using a linear interpolation
c between 2 regions of constant Cv
c
c John Mahaffy 2/3/96
c
implicit none
real Cv,u,T,w
c
c Cv - specific heat at constant volume
c u - internal energy per kilogram
c T - temperature
c w - interpolation weighting function
c
T=300.
do 100 while (T.le.3000.)
if(T.lt.1600.) then
Cv=5000.
else if (T.gt.2200.) then
Cv=7000.
else
w=(T-1600.)/(2200.-1600.)
Cv= (1.-w)*5000.+w*7000.
endif
u=Cv*T
write(10,*)T,Cv,u
T=T+10.
100 continue
return
end
subroutine cubenrgy
c
c Generate internal energy and Cv of a gas for Temperatures
c between 300 and 3000K in steps of 10K using a cubic interpolation
c between 2 regions of constant Cv.
c The cubic is chosen so that derivatives are continuous.
c
c John Mahaffy 2/3/96
c
implicit none
real Cv,u,T,ST,w
c
c Cv - specific heat at constant volume
c u - internal energy per kilogram
c T - temperature
c ST - Scaled temperature
c w - interpolation weighting function
c
T=300.
do 200 while (T.le.3000.)
if(T.lt.1600.) then
Cv=5000.
else if (T.gt.2200.) then
Cv=7000.
else
ST=(T-1600.)/(2200.-1600.)
w= ST**2*(3.-2.*ST)
Cv= (1.-w)*5000.+w*7000.
endif
u=Cv*T
write(11,2000)T,Cv,u
T=T+10.
200 continue
2000 format(1x,f7.1,1p,2e12.5)
return
end
subroutine plotit
c
c This subroutine is for those of you who are getting bored. At this
c stage of the class you will not be held responsible new material here.
c
c<a name="logical"><font color="FF0000">
logical fexist
c</font></a>
c fexist - flag indicating whether or not the gnuin file exists
c First check to see if a file containing commands to the graphic program
c exists in your local directory
c
inquire (file='gp-int1',exist=fexist)
c
c If it doesn't then issue a unix command to copy it from my directory
c the subroutine "system" is not a standard part of Fortran, but is
c provided by some Unix systems to give a crude form of what is called
c a "system call". Any computer will supply an equivalent command, and
c other functions or subroutines that give more direct access to system
c support.
c
if (.not.fexist) then
call system('cp ~jhm/201/gp-int1 .')
print *, ' Graphics input file copied'
endif
c
c Close the data files so they are ready for reading by gnuplot
c
close ( 11 )
close ( 10 )
c
c If you are sitting at a Hammond lab work station remove the "c" in
c column 1 of the following "call system" to plot your results. If you
c are running from NCSA/BYU Telnet on a Mac, you can see plots by
c making the same change, and by following instructions in "gp-int1"
c to change that file for Tektronics emulation in gnuplot.
c
c call system ('gnuplot gp-int1')
c
return
end
c</pre>
c</body>
c</html>