program heatflow; {by Steve Smith. Simulates heat flow through concentric
water or glass rings of thickness dr and width dL. Heat power of 'input'
watts/meter is deposited according to Beers law inside cylinder of radius
'rBeam'. This version implements lateral heat flow and removes radial
temperature cutoff. Temperature data are saved in OUTPUT folder as a text
file and as two bitmap images.}
uses crt, graph, dos;
const
fname='MYNAME'; {change filename for new data output file, 6 letters max}
dt=4e-6; {iteration time in seconds}
input= 19.5e-3;{Watts  =   111 mW into objective lens transmitting 17.6%}
absCoeff = 3140; {absorbance coeff. of water per meter.....use natural logs}
dr=5e-6; {radial step size in meters}
dL=5e-6; {lateral step size in meters}
iMax = 70; {70 x 70 maximum array size}
jMax = 70; {array goes from -jMax to +jMax so it counts twice as large}
rBeam=30e-6; {radius of heating beam in meters. Use multiples of dr}
innerWall= 55e-6; { +/- inside lateral distance to chamber walls }
outerWall= 255e-6;{+/- distance to outside chamber, assumes #2 coverglass}
waterCond = 0.6 ; {W/deg-meter}   {heat conductivity of water}
glassCond = 1.1 ; {W/deg-meter}   {borosilicate glass}
CargCond  = 0.1 ; {W/deg-meter}   {Cargille seriesAAA n=1.33}
waterSpHeat = 4.18e6; {J/meter^3-deg}  {volume-specific heat capacity}
glassSpHeat = 1.67e6; {J/meter^3-deg}
CargSpHeat =  1.61e6; {J/meter^3-deg}

type bigArray = array[0..iMax+1, -jMax..jMax] of real;
{stores present and future temperatures at cylindrical coordinates r and L}

var
tempNow: bigArray;  {will reside in data segment}
tempNext:^bigArray; {will reside in heap}
i,j,k,xMax,yMax,xMid,yMid,x,y,color:integer;
time,fluxIn,fluxOut,deltaFlux,r,L:real;
ch:char; outFile:text;

{$I BITMAP.PAS}

procedure bigspot(x,y,color:integer);
var i,j:integer;
begin
for i:=-2 to 2 do
for j:=-2 to 2 do
putPixel(x+i,y+j,color);
end;

procedure showMap; {temperature map of plane thru center of Kohler beam}
var i,j,color,idR,jdL:integer;
begin
for i:=1 to iMax do for j:= -jMax to jMax do
  begin    {color coded temperature map in r-L coordinates}
  color:=round(tempNow[i,j]);
  idR:=round(i*dR*1e6); jdL:= round(j*dL*1e6);
  bigspot(xMid+idR-1,yMid-jdL+1,color);
  bigspot(xMid-idR+1,yMid-jdL+1,color); {duplicate symmetric halves}
  end;
setColor(white);      {white lines indicate glass chamber walls}
line(xMax,yMid-round(innerWall*1e6),0,yMid-round(innerWall*1e6));
outTextXY(50,100,'GLASS');
line(xMax,yMid+round(innerWall*1e6),0,yMid+round(innerWall*1e6));
outTextXY(50,yMax-100,'GLASS');
outTextXY(50,yMid,'WATER');
setTextStyle(0,1,1);
outTextXY(xMid,0,'heating laser beam');
setTextStyle(0,0,1);
setColor(red);        {red outlines heating-laser beam}
line(xMid+round(rBeam*1e6),yMax,xMid+round(rBeam*1e6),0);
line(xMid-round(rBeam*1e6),yMax,xMid-round(rBeam*1e6),0);
ch:=' ';
end;

procedure saveEverything;
var i,j:integer; message:string;
begin
assign(outFile,'output\'+fname+'.TXT');
rewrite(outFile);
writeln(outFile,'HEATFLOW.PAS      dt = ', dt:3,' seconds');
writeln(outFile,'Input = ',input*1000:1:2,' mW     absCoeff = ',absCoeff);
writeln('absCoeff=',absCoeff);
writeln(outFile,'rBeam = ',rBeam*1e6:1:1,' microns');
writeln(outFile,'innerWall = +/-',innerWall*1e6:1:1,' microns');
writeln(outFile,'time = ',time*1000:1:2,' ms');
for i:= 0 to iMax+1 do write(outFile,chr(9),'r=',round(i*dr*1e6));
writeln(outFile);
for j:=-jMax to jMax do
  begin
  write(outFile,'L=',round(j*dL*1e6),chr(9));
  for i:=0 to iMax do write(outFile,tempNow[i,j]:1:3, chr(9));
  writeln(outFile);
  end;
close(outFile);
assign(outFile,'output\'+fname+'.BMP');
rewrite(outFile);
writeBitmap;   {save bitmap image of temperature vs R, and risetime history}
close(outFile);
clearDevice;
showMap;
assign(outFile,'output\'+fname+'2.BMP');
rewrite(outFile);
writeBitmap;
close(outFile);
ch:='q';
end;

procedure showMenu(hue:word);
var message:string;
begin
setColor(hue);
str(tempNow[0,0]:1:3,message);
outTextXY(xMid,yMid-40,'temp[0,0,now] = '+message+' deg');
str(time*1000:1:1,message);
outTextXY(xMid,yMid-20,'time = '+ message+' ms');
outTextXY(xMid,yMid,   'press spacebar to continue');
outTextXY(xMid,yMid+20,'s = save         q = quit');
end;

function conductanceR(i,j:integer):real;
{Radial conductance of ring[i,j], units are Watts/deg}
var r,L,cond:real;
begin
r:=dr*i;  L:=dL*j;
if abs(L)>outerWall then cond:=CargCond
   else if abs(L)>innerWall then cond:= glassCond
        else cond:= waterCond;
conductanceR := cond*2*Pi*r*dL / dr;
end;

function conductanceL(i,j:integer):real;
{Lateral conductance of ring[i,j], units are Watts/deg}
var r,L,cond:real;
begin
r:=dr*i;  L:=dL*j;
if abs(L)>outerWall then cond:=CargCond
   else if abs(L)>innerWall then cond:= glassCond
        else cond:= waterCond;
conductanceL:= cond*2*Pi*r*dr / dL;
end;

function capacitance(i,j:integer):real;
{Amount of energy to raise temperature of ring[i,j] by 1 degree
Output units are Joules/deg}
var r,L,spHeat:real;
begin
r:=dr*i;  L:=dL*j;
if abs(L)>outerWall then spHeat:=CargSpHeat
   else if abs(L)>innerWall then spHeat:= glassSpHeat
        else spHeat:= waterSpHeat;
capacitance := spHeat*Pi*dL*(2*r*dr-dr*dr);
end;

function heat(j:integer):real; {units of power, same as INPUT}
var L,lightPower:real;
begin
L:= j*dL;
if abs(L)<innerWall { L runs both negative and positive }
then begin              {remove negative offset for Beers law}
     lightPower:=input*exp(-absCoeff*(L+innerWall));
     heat:=lightPower*absCoeff*dL;
     end
else heat:=0;
end;

begin  {main}
New(tempNext);
i:=detect;
initgraph(i,j,'');
xMax:=getMaxX; yMax:= getMaxY; xMid:=xMax div 2; yMid:=yMax div 2;
repeat
for j:=-jMax to jMax do
    begin
    for i:=1 to iMax do
        begin
        r:=dr*i;  L:=dL*j;
        fluxIn :=  (tempNow[i-1,j]-tempNow[i,j])
                 * (conductanceR(i-1,j)+conductanceR(i,j))/2
                 + (tempNow[i,j-1]-tempNow[i,j])
                 * (conductanceL(i,j-1)+conductanceL(i,j))/2;
        fluxOut:=  (tempNow[i,j]-tempNow[i+1,j])
                 * (conductanceR(i+1,j)+conductanceR(i,j))/2
                 + (tempNow[i,j]-tempNow[i,j+1])
                 * (conductanceL(i,j)+conductanceL(i,j+1))/2;
        deltaFlux := fluxIn-fluxOut;
        if (r<=rBeam) then deltaFlux:=deltaFlux+Heat(j)*2*r*dr/(rBeam*rBeam);
        tempNext^[i,j] := tempNow[i,j] + deltaFlux * dt / capacitance(i,j);
        tempNext^[i, jMax] :=0; {lateral cutoffs}
        tempNext^[i,-jMax] :=0; {at both ends, which are objective lenses}
        end; {for i}
    tempNext^[0,j]:=tempNext^[1,j];    {reflecting wall for mirror symmetry}
    end; {for j}
for i:=0 to iMax+1 do
for j:=-jMax to jMax do
    begin
    tempNow[i,j]:=tempNext^[i,j];  {update temperatures}
    x:=i*round(200/iMax) ;
    y:= yMax-round(tempNow[i,j]*10);
    color:= round(tempNow[i,j]);
    if k mod 10 = 0 then
       begin
       putPixel(x+0,y,color);        {bar graph of temperature vs radius}
       putPixel(x+1,y,color);
       putPixel(x+2,y,color);
       putPixel(x+3,y,color);
       end;
    end;
time:=time+dt;
inc(k);
putPixel(k div 25 +200, yMax-round(tempNow[0,0]*10),cyan);
if (k>2500) then            {temperature-rise vs time plot}
   begin
   showMenu(white);
   ch:=readkey;
   k:=0;
   end;
if ch = 's' then saveEverything;
if k=0 then showMenu(black);
until (ch='q');
dispose(tempNext);
end.