c
c hotspot - model heat dissipation by conduction and radiation
c           in two dimensions under a localized heat load
c
c Richard T. Jones, July 23 2000
c--------------------------------

c Model
c
c The model consists of a thin film characterized by 
c    1) the film thickness (length)
c    1) heat capacity Cp (Joules/K/volume)
c    2) thermal conductivity kappa (Watts/K/length)
c    3) Steffan-Boltzmann radiation constant SBc (Watts/area/K**4)
c The driving term in the model is the heating function h(x,y) that is
c expressed in Watts/area.  The basic unit for length is the cm, sorry
c for mixing cgs and mks units.  Time in this model is artificial - only
c the steady state is of interest - so the value of Cp is adjusted
c to give a gradual time evolution using unit time steps.  All of the
c other constants are physical.  The equation of motion is as follows.
c
c      dT         1                      4   4                  _2
c      --  =  --------- ( h(x,y) - SBc (T - T )  + kappa rthick V T )
c      dt     rthick Cp                      0
c
c The funny V with a _2 overhead is the LaPlacian.  Boundary conditions
c are that the temperature is held to Tambient at the perimeter of the
c grid, which is a rectangular shape.

c Procedure
c
c 1) call initialize() to initialize the thermal grid.  This will set
c    the temperature everywhere to a uniform value given by Tambient.
c 2) call hotspot(power) to set the total power level of the heater.  This
c    can be called at any time to change the heating term.
c 3) call diffuse(nsteps) to advance nsteps time units forward.  The state
c    of the thermal grid can be measured at any time.



      subroutine hotspot(current)
      implicit none
      real current
      include 'hotspot.inc'
      data rthick,rwide,rhigh /15.e-4,0.5,0.5/
      data beamsizeX,beamsizeY,hpower/0.17,0.05,10.e-3/
      real BetheBloch
      external BetheBloch
      real area
      area = 2*pi_*beamsizeX*beamsizeY
      hpower = current*BetheBloch(12e9)*density*rthick
      print *, 'Power on crystal in Watts is',hpower
      hpower = hpower/area
      end


      subroutine initialize()
      implicit none
      include 'hotspot.inc'
      integer ix,iy
      time = 0
      do ix=1,NX
        do iy=1,NY
          Tgrid(ix,iy) = Tambient
        enddo
      enddo
      end


      subroutine diffuse(nsteps)
      implicit none
      integer nsteps
      include 'hotspot.inc'
      integer i
      do i=1,nsteps
        call diffuse2
      enddo
      end


      subroutine diffuse1        ! with thermal clamps at periphery
      implicit none
      include 'hotspot.inc'
      real T(NX,NY),dT
      real x,y
      integer ix,iy
      do ix=1,NX
        do iy=1,NY
          T(ix,iy) = Tgrid(ix,iy)
        enddo
      enddo
      time = time+1
      do ix=2,NX-1
        x = ((ix-0.5)/NX-0.5)*rwide
        do iy=2,NY-1
          y = ((iy-0.5)/NY-0.5)*rhigh
          dT = 1./(Cp*rthick) * (
     +           hpower * exp(-0.5*(
     +             (x/beamsizeX)**2 + (y/beamsizeY)**2)
     +           )
     +         - 2*SBc*(T(ix,iy)**4 - Tambient**4)
     +         + kappa*rthick * (
     +             (T(ix-1,iy) + T(ix+1,iy) - 2*T(ix,iy))
     +               /(rwide/NX)**2
     +           + (T(ix,iy-1) + T(ix,iy+1) - 2*T(ix,iy))
     +               /(rhigh/NY)**2
     +           )
     +         )
          Tgrid(ix,iy) = Tgrid(ix,iy) + dT
        enddo
      enddo
      end


      subroutine diffuse2        ! insulated, only radiative cooling
      implicit none
      include 'hotspot.inc'
      real T(NX,NY),dT
      real x,y
      integer ix,iy
      vector dTmax(1)
      do ix=1,NX
        do iy=1,NY
          T(ix,iy) = Tgrid(ix,iy)
        enddo
      enddo
      time = time+1
      dTmax(1) = 0
      do ix=2,NX-1
        x = ((ix-0.5)/NX-0.5)*rwide
        do iy=2,NY-1
          y = ((iy-0.5)/NY-0.5)*rhigh
          dT = 1./(Cp*rthick) * (
     +           hpower * exp(-0.5*(
     +             (x/beamsizeX)**2 + (y/beamsizeY)**2)
     +           )
     +         - 2*SBc*(T(ix,iy)**4 - Tambient**4)
     +         + kappa*rthick * (
     +             (T(ix-1,iy) + T(ix+1,iy) - 2*T(ix,iy))
     +               /(rwide/NX)**2
     +           + (T(ix,iy-1) + T(ix,iy+1) - 2*T(ix,iy))
     +               /(rhigh/NY)**2
     +           )
     +         )
          Tgrid(ix,iy) = Tgrid(ix,iy) + dT
          dTmax(1) = max(dTmax(1),dT)
        enddo
        Tgrid(ix,1) = Tgrid(ix,2)
        Tgrid(ix,NY) = Tgrid(ix,NY-1)
      enddo
      do iy=1,NY
        Tgrid(1,iy) = Tgrid(2,iy)
        Tgrid(NX,iy) = Tgrid(NX-1,iy)
      enddo
      end


      subroutine diffuse3        ! insulated along vertical boundaries
      implicit none
      include 'hotspot.inc'
      real T(NX,NY),dT
      real x,y
      integer ix,iy
      do ix=1,NX
        do iy=1,NY
          T(ix,iy) = Tgrid(ix,iy)
        enddo
      enddo
      time = time+1
      do ix=2,NX-1
        x = ((ix-0.5)/NX-0.5)*rwide
        do iy=2,NY-1
          y = ((iy-0.5)/NY-0.5)*rhigh
          dT = 1./(Cp*rthick) * (
     +           hpower * exp(-0.5*(
     +             (x/beamsizeX)**2 + (y/beamsizeY)**2)
     +           )
     +         - 2*SBc*(T(ix,iy)**4 - Tambient**4)
     +         + kappa*rthick * (
     +             (T(ix-1,iy) + T(ix+1,iy) - 2*T(ix,iy))
     +               /(rwide/NX)**2
     +           + (T(ix,iy-1) + T(ix,iy+1) - 2*T(ix,iy))
     +               /(rhigh/NY)**2
     +           )
     +         )
          Tgrid(ix,iy) = Tgrid(ix,iy) + dT
        enddo
      enddo
      do iy=1,NY
        Tgrid(1,iy) = Tgrid(2,iy)
        Tgrid(NX,iy) = Tgrid(NX-1,iy)
      enddo
      end


      real function BetheBloch(E)        ! in eV
      implicit none
      real E
      include 'hotspot.inc'
      real Tmax,Tcut,beta,gamma,del
      real hbwp
      Tmax = E/2
      Tcut = min(Tmax,15.e6*density*rthick)
      gamma = 1+E/me_c2
      beta = sqrt(1-1/gamma**2)
      hbwp = 28.816*sqrt(density*Z_/A_)
      del = 2*(log(hbwp/I_)+log(beta*gamma)-0.5)
      BetheBloch = K_*(Z_/A_)/beta**2 * (
     +               0.5*log(2*me_c2*Tcut*(beta*gamma/I_)**2)
     +             - (beta**2)/2*(1+Tcut/Tmax)
     +             - del/2
     +             )
      end