c c c
      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
       logical fexist
c

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
c c