#!/usr/bin/perl
require "code.gplo";

make_f77(); system("./a.out  ");
#exit;

#....... make plots, loop over output files
$k=0;
foreach $name (glob("tmpname-temp*")){  if($name=~/(\d{6})/){ $ist=$1; 
   $pl1=" xlab=depth(m) lwidth=3 type=1 with=lines ";   $ou='tmpx.pdf';
   $ylo='set~logscale~y;set~format~y~"10^{%L}"';
   $gpl=gpl_head("term=jpg font=Helvetica fsize=6 dx=0.5 title=1d~heat~equation~$ist");
   $gpl=$gpl . gpl_xyplot("file=$name using=1:2 key=Temperature(C) yfix=0:100  ylab=Temperature $pl1 ");
   #$gpl=$gpl . gpl_xyplot("file=$name using=1:3 key=Conductivity ylab=kappa(W/m/K) lcol=rgb~'#dd9999' $pl1  ");
   #$gpl=$gpl . gpl_xyplot("file=$name using=1:4 key=Source ylab=Power(W/m^3) lcol=rgb~'#5555ff' yfix=1.0e-9:0.9 $pl1 more=$ylo ");
   gpl_final($gpl,"out=$ou");
   ssy("convert -quality 100 -crop 800x400+0+540  tmp.jpg tmp_$k.jpg");  $k++;
   #  ssy("xv tmp_$ist.jpg");
}}
#....... make movie
system("/bin/rm -fv movie.mp4");  #...... clean out stale files
system("ffmpeg  -i tmp_%d.jpg  movie.mp4");
print STDERR "\n\n  mpv --loop=inf movie.mp4\n\n";

sub ssy{ local($c)=@_; print STDERR "ssy:|$c|\n"; system($c); return; }

sub make_f77{
$t="
c---------------------------------------------------
      program heat1d
c---------------------------------------------------
c     time-dependent heat equation
c     JR  11/2017

c     heat equation for vertical earth profile
c     x(i) is depth in meter
c     u(i) is temperature in C

      external xkappa,rho
      parameter(N=900)
      real x(N),T(N),Tnew(N),d(N),a(N),b(N),c(N)
      character*80 fname

      time0=0.0
      yy=3600.0*24.0*365.0
      timef=1000.0*yy
      x1=0.0
      t1=10.0     !...... left boundary, earth surface
      xn=3000.0
      tn=19.0    !...... right boundary, 3km depth, assuming 3k/km gradient
      dx=(xn-x1)/float(N-1)
      dt=7.5    !.... unstable
      dt=5.5    !.... stable
      dt=5.6    !.... unstable
      cfl=0.5*dx*dx
      ck=xkappa(1000.0)/rhocv(1000.0)
      cfl=cfl/ck    !..... adjust for real values
      write(0,*)'cfl,dt  ',cfl/yy,dt
      dt=0.1*yy
      dtout=10.0*yy
       ! stop
c...... initial values, step function
      do i=1,N
          x(i)=x1+float(i-1)*dx
          !  T(i)=t1
          !  if(i.gt.N/2)T(i)=tn
          T(i)=t1+float(i-1)*(tn-t1)/float(N-1)
          Tnew(i)=T(i)
      enddo
c...... time step look
      time=time0
      kstep=0
      outnext=time0
100   continue
c..... update interior interior points, Euler
      do i=2,N-1
          xx=x(i)
          xm=xx-0.5*dx
          xp=xx+0.5*dx
         Tnew(i)=T(i)
     +           + (dt/dx/dx/rhocv(xx))
     +             * (  xkappa(xm)*T(i-1)
     +                + xkappa(xp)*T(i+1)
     +                - (xkappa(xm)+xkappa(xp))*T(i) )
     +                + dt*qsource(x(i),time)
       !   write(0,*)i,x(i),rhocv(xx),xkappa(xm),qsource(x(i),time)
       !   write(0,*)i,T(i),Tnew(i)
       if(Tnew(i).gt.1000000000.0) then
         write(0,*)'solution is blowing up ',i
         stop
       endif
      enddo
      !...... copy back
      do i=1,n
        T(i)=Tnew(i)
      enddo

      time=time+dt
      !........  output some points to look at progress
      write(0,*)time/yy,T(300),T(440),T(455),T(600)

      !...... output every so often
      if(time.gt.outnext) then
        k=int(dtout+0.5); m=int(time); m=m/k; iout=m*k;
        write(0,*)' m,k,iout ',m,k,iout
        fname=' '
        write(fname,'(a,i6.6)')'tmpname-temp.',int(time/yy) 
        call wout(n,x,T,fname)
        outnext=time+dtout
      endif
      !...... next step if still time left
      kstep=kstep+1
      if(kstep.gt.10000)stop
      if(time.lt.timef) goto 100

      stop
      end
c---------------------------------------------------
      subroutine wout(n,x,y,fname)
c---------------------------------------------------
      real x(n),y(n)
      character*(*) fname
      write(0,*)'writing ... ',TRIM(fname)
      open(10,file=fname)
      do i=1,n
      write(10,*)x(i),y(i)
      enddo
      close(10)
      return
      end
c---------------------------------------------------
      subroutine thomas(n,a,b,c,d,x)
c---------------------------------------------------
      real a(n),b(n),c(n),d(n),x(n)  !... note we overwrite c and d
      c(1)=c(1)/b(1)
      d(1)=d(1)/b(1)
      do i=2,n
          z=1.0/(b(i)-a(i)*c(i-1))
          c(i)=c(i)*z
          d(i)=(d(i)-a(i)*d(i-1))*z
      enddo
      x(n)=d(n)
      do i=n-1,1,-1
          x(i)=d(i)-c(i)*x(i+1)
      enddo
      return
      end
c---------------------------------------------------
      real function rhocv(x)  
c---------------------------------------------------
c   returns material rho*Cv
      cv=2000    !..... heat capacity in J/K/kg
      rho=2750   !..... density of granite in kg/m^3
      rhocv=rho*cv   !.....  in J/K/m^3
          !rhocv=1.0  !...... for testing and stbility demo
      return
      end
c---------------------------------------------------
      real function xkappa(x) 
c---------------------------------------------------
c   returns thermal conductivity in W/m/K
      xkappa=1.8                       !... granite
      if(x.lt.1000.0) xkappa=2.2       !.... sandstone
      if(x.lt.500.0)  xkappa=1.0       !.... limestone
      if(x.lt.150.0)  xkappa=0.5       !.... soil/gravel 
!      xkappa=1.0                       !.... testing and demo
      return
      end
c---------------------------------------------------
      real function qsource(x,time)
c---------------------------------------------------
c.....   returns heating rate in W/m^3
      yy=3600.0*24.0*365.0
       qsource=0.0      
      !    qsource=1.0e-8
      ! if(x.gt.1000.0)qsource=0.5e-6                     !...... granite
      if((x.gt.800.0).and.(x.lt.825.0)) then
        qsource=1.0e-7*exp(-time/yy/30.0)     !...... diluted radioactive waste after 10s of years
        write(0,*) 'qsource: ',x,qsource
      endif
!      qsource=0.0
      return
      end
    ";

    system("/bin/rm -f tmp.f a.out");          #...... clean up
    open(FF,">tmp.f"); print FF $t; close(FF);    #...... write code to temporary file
    system("gfortran tmp.f");
    if( -e "a.out" ){ system("  size a.out  "); return; }
    else{ print STDERR "compilation failed!\n"; exit; }
}