(* 
Package for calculation in cryogenics for Mathematica
<C> Rok Zitko 2004, rok.zitko@ijs.si, maj 2004
$Id: krio.m,v 1.3 2005/02/04 10:19:56 rok Exp rok $
*)

(*
CHANGE LOG

2.2.2005 - system simulations, quasistationary systems only!
3.2.2005 - AmountBy* - units changed!!
4.2.2005 - ShieldSolve, digitized Conductivity data
5.2.2005 - RoundPlateSolve, usage messages

*)

(* 
TO DO

Thermal sources/sinks
*)

(* Turn off annoying warnings *)
Off[General::spell1];
Off[General::spell];

BeginPackage["Krio`"];

Needs["LinearAlgebra`MatrixManipulation`"];
Needs["Statistics`DataManipulation`"];

(* Convenience definitions *)
m=1; cm=m/100; mm=m/1000; um=mm/1000;
l=1/1000; (* liter *)

(* 
NOTE: All values are expressed in standard SI units throughout
the package! 
*)

(* Source: cryo_materials.pdf, digitized using Engauge Digitizer *)

ConductivityList::"usage"="ConductivityList[m] returns tabulated
data for temperature dependant thermal conductivity for material m.";

ConductivityList[Cu100]=
{
{1.01418,144.856},
{1.97564,304.449},
{3.04135,482.199},
{3.95042,651.305},
{6.00231,1013.49},
{8.00193,1345.07},
{9.99268,1663.2},
{14.9819,2328.02},
{19.4503,2635.83},
{24.9121,2501.51},
{29.89,2251.27},
{34.9224,1698.22},
{39.2531,1475.1},
{50.2279,1001.18},
{61.0245,755.308},
{69.4851,633.344},
{82.2612,503.744},
{93.7035,486.439},
{102.65,478.043},
{121.578,445.66},
{138.503,445.808},
{157.769,430.494},
{186.899,430.68},
{245.752,438.653},
{294.946,438.858}
};

ConductivityList[Cu20]=
{
{1.01491,25.4456},
{1.95027,56.2341},
{3.00107,90.4967},
{3.84704,131.023},
{5.84337,214.601},
{7.5868,269.85},
{9.85171,364.103},
{14.7688,565.654},
{19.4274,723.926},
{24.5715,848.343},
{29.494,910.298},
{34.0418,926.481},
{38.7756,863.425},
{49.0168,763.226},
{58.8112,651.292},
{66.1241,617.755},
{76.2997,546.066},
{89.1933,474.264},
{100.284,449.843},
{120.351,434.265},
{137.105,434.265},
{160.305,419.227},
{189.903,419.227},
{249.698,426.68},
{299.672,419.227}
};

ConductivityList[Al]={
{1.02607,8.38435},
{1.99712,16.3789},
{3.07286,25.0011},
{3.93881,34.9436},
{6.06084,55.2519},
{7.97265,70.7115},
{9.95477,90.4967},
{14.9223,135.723},
{19.63,176.786},
{24.8285,210.853},
{29.4196,238.534},
{34.8573,260.505},
{39.2021,284.5},
{49.5691,289.557},
{58.7194,284.5},
{69.5449,251.484},
{80.2656,251.484},
{89.0807,242.775},
{100.167,242.775},
{120.203,226.252},
{136.927,218.417},
{155.989,218.417},
{175.391,210.853},
{197.22,210.853},
{249.366,210.853},
{303.21,210.853}
};

ConductivityList[SS]={
{1.01578,0.0402339},
{2.00397,0.100589},
{3.04449,0.183128},
{3.9539,0.265136},
{6.00607,0.449843},
{7.90245,0.651292},
{9.99756,0.878775},
{14.9884,1.4142},
{19.7209,2.0475},
{24.9485,2.7144},
{29.9556,3.41321},
{35.0386,4.07093},
{39.4073,4.52493},
{49.8484,5.68987},
{60.6279,6.55129},
{69.9926,7.54312},
{79.7548,8.53341},
{92.0586,8.99667},
{102.186,9.48508},
{121.073,10.3587},
{139.756,11.1152},
{159.232,11.927},
{181.405,12.139},
{196.174,12.5745},
{254.628,13.4928},
{297.775,14.4782}
};

(* Chapter 3: Cryogenic heat transfer *)

(* Thermal conductivity integrals *)
(* Form: {{T_1 [K],int K dT [W/m]},...} *)

ConductivityIntegralList::"usage"="ConductivityIntegralList[m] returns
tabulated thermal conductivity integral for material m.";

ConductivityIntegralList[Cu100]=
{
{300,1.8 10^5},{200,1.5 10^5},{100,1.07 10^5},
{70,9.7 10^4},{60,8.8 10^4},{50,7.9 10^4},
{40,6.7 10^4},{30,5.1 10^4},{20,2.9 10^4},
{10,7.6 10^3},{5,1.8 10^3},{1,7.1 10}
};
ConductivityIntegralList[Cu20]=
{
{300,1.34 10^5},{200,9.2 10^4},{100,5.8 10^4},
{90,5.4 10^4},
{70,4.4 10^4},{50,3 10^4},{30,1.5 10^4},
{20,6.8 10^3},{10,1.6 10^3},{4,2.25 10^2},
{2,5.1 10^1},{1,1.2 10^1}
};
ConductivityIntegralList[SS]=
{
{300,2.9 10^3},{200,1.5 10^3},{100,5.1 10^2},
{80,3.3 10^2},{50,1.3 10^2},{40,8.1 10^1},{30,4.4 10^1},{20,1.7 10^1},
{10,3.8},{5,7 10 ^-2},{1,3.9 10^-2}
};
ConductivityIntegralList[PTFE]=
{
{300,4.7 10^1},{200,2.9 10^1},{100,1.3 10^1},
{80,1.07 10^1},{60,5.8 10^0},{40,3.6 10^0},
{20,1.4 10^0},{15,1.0 10^0},{10,3.4 10^-1},
{8,2.25 10^-1},{6,1.2 10^-1},{4,3.4 10^-2},
{2,5.1 10^-3},{1,8.7 10^-4}
};

(* Load time interpolations *)
SetOptions[Interpolation,InterpolationOrder->1];
ConductivityIntegral[Cu100]=Interpolation[ConductivityIntegralList[Cu100]];
ConductivityIntegral[Cu20]=Interpolation[ConductivityIntegralList[Cu20]];
ConductivityIntegral[SS]=Interpolation[ConductivityIntegralList[SS]];
ConductivityIntegral[PTFE]=Interpolation[ConductivityIntegralList[PTFE]];
Conductivity[Cu100]=Interpolation[ConductivityList[Cu100]];
Conductivity[Cu20]=Interpolation[ConductivityList[Cu20]];
Conductivity[SS]=Interpolation[ConductivityList[SS]];
Conductivity[Al]=Interpolation[ConductivityList[Al]];

(* Run time macros for conduction integrals *)

ConductivityIntegral::"usage"="ConductivityIntegral[material, T] returns
the value of conductivity integral for material at temperature T.
ConductivityIntegral[material, {T1, T2}] returns the conductivity
integral between temperature T1 and T2, as needed for calculation of 
heat conductivity";

ConductivityIntegral[material_,T_?NumberQ]:=ConductivityIntegral[material][T];
ConductivityIntegral[material_,T_List /; Length[T]==2]:=
ConductivityIntegral[material][T[[1]]]-ConductivityIntegral[material][T[[2]]];

(* Emissivity of surfaces *)
Emissivity::"usage"="Emissivity[material] returns the emissivity of
material. Note that it's hard to perform measurements of emissivity, so these
values should only be considered as a rough guide.";

Emissivity[CuPolished]=0.024;
Emissivity[CuOxydized]=0.78;
Emissivity[AgPolished]=0.02;
Emissivity[AuPolished]=0.018;
Emissivity[SSPolished]=0.11;
Emissivity[SSSanded]=0.38;

(* Heat conduction formulas *)
ConductionHeatTransfer::"usage"="ConductionHeatTransfer[G, m, {T1, T2}]
returns the heat current through an object made of material m 
between temperatures T1 and T2. The geometry of the object is 
described using the geometry factor G. See function GeomRod,
GeomPipe, GeomSquareRod, GeomBaffle.";

ConductionHeatTransfer[geomfactor_,material_,tmps_]:=
geomfactor ConductivityIntegral[material,tmps];

(* Heat radiation formulas *)
StefanConstant=5.7 10^-8; (* W/(m^2K^4) *)
HeatRadiationCylinders::"usage"="HeatConductionCylinders[A1, A2, eps1, eps2, T1, T2]
returns the radiative heat current between long concentric cylindrical 
surfaces with surface areas A1 and A2. The emissivities are given by
eps1 and eps2 and the temperatures by T1 and T2. 1 refers to the
inner (colder?) surface, 2 refers to the outer (hotter?) surface.
Note that maximum radiative heat current from a surface at room
surface to a surface at LHe (???) temperature is 46mW/cm^2, while
the maximum current from a surface at LN2 temperature to a surface 
at LHe temperature is 300microW/cm^2";
HeatRadiationCylinders[A1_,A2_,eps1_,eps2_,T1_,T2_]=
A1 StefanConstant (T2^4-T1^4)/(1/eps1+A1/A2(1/eps2-1));

HeatRadiationPlane::"usage"="HeatRadiationPlane[a, eps1, eps2, T1, T2]
returns the radiative heat current between two plane surfaces of area a
with emissivities eps1 and eps2 and temperatures T1 and T2.";
HeatRadiationPlane[area_, epsilon1_, epsilon2_, T1_, T2_]=
StefanConstant area epsilon1 epsilon2/(epsilon1+epsilon2-epsilon1 epsilon2) (T2^4-T1^4);
(* White p. 220 *)

(* Cross sections *)
CrossSectionPipe::"usage"="CrossSectionPipe[od, id] returns area of
cross section of a circular hollow pipe with outer diameter od and 
inner diameter id.";
CrossSectionPipe[od_,id_]:=Pi ((od/2)^2-(id/2)^2);

CrossSectionRod::"usage"="CrossSectionRod[r] returns cross section
area of a rod of radius r.";
CrossSectionRod[r_]:=Pi r^2;

CrossSectionSquareRod::"usage"="CrossSectionSquareRod[a] returns cross
section area of a square rod of side a.";
CrossSectionSquareRod[a_]:=a^2;

(* Surface Areas *)
SurfaceAreaCylinder::"usage"="SurfaceAreaCylinder[r, l] returns surface
area of cylinder of radius r and length l.";
SurfaceAreaCylinder[r_,l_]:=2 Pi r l;

SurfaceAreaTube::"usage"="SurfaceAreaTube[r, l] returns surface area
of a cylindical tube of radius r and length l that is closed on one
end.";
SurfaceAreaTube[r_,l_]:=2 Pi r l + Pi r^2;

SurfaceAreaTubeClosed::"usage"="SurfaceAreaTube[r, l] returns surface
area of a cylindrical tube of radius r and length l that is closed on
both ends.";
SurfaceAreaTubeClosed[r_,l_]:=2 Pi r l+2 Pi r^2;

(* Geometry factor G, 1/Integrate[dx/A(x)] *)
GeomPipe::"usage"="GeomPipe[od, id, l] returns the geometry factor
of a hollow circular pipe with outer diameter od, inner diameter id 
and length l.";
GeomPipe[od_,id_,length_]:=CrossSectionPipe[od,id]/length;

GeomRod::"usage"="GeomRod[r, l] returns the geometry factor
of a round rod of radius r and length l.";
GeomRod[radius_,length_]:=CrossSectionRod[radius]/length;

GeomSquareRod::"usage"="GeomSquareRod[a, l] returns the geometry
factor for a square rod of side a and length l.";
GeomSquareRod[a_,l_]:=CrossSectionSquareRod[a]/length;

GeomBaffle::"usage"="GeomBaffle[r2, r1, t] returns the
geometry factor for radiation baffle of outer radius r2,
inner radius r1 and thickness t.";
GeomBaffle[r2_, r1_, thickness_]:=2 Pi thickness/Log[r2/r1];

(* Constants *)
GasConstant=8.31441; (* J/(mol K) *)

AmountBy1W::"usage"="AmountBy1W[l] returns the volume of
cryogenic liquid l evaporated in one hour by 1 Watt of power.";
AmountBy1W[N]=0.024/1000; (* m^3/hour, for 1 Watt of power *)
AmountBy1W[He]=1.38/1000; (* m^3/hour *)

(* Common temperatures *)
TempRT=293;
TempLN=77.4;
TempLHe=4.2;
RTLN={TempRT,TempLN};
RTLHe={TempRT,TempLHe};
LNLHe={TempLN,TempLHe};

(* Liquid densities *)
LiquidDensityNBP::"usage"="LiquidDensityNBP[l] returns the
density of liquid l at its normal boiling point.";
LiquidDensityNBP[N]=808; (* kg/m^3 *)
LiquidDensityNBP[He]=125;

GasDensityNBP::"usage"="GasDensityNBP[g] returns the 
density of gas l at its normal boiling point.";
GasDensityNBP[N]=1.16;
GasDensityNBP[He]=0.166;

(* Normal boiling points, Weisend p. 7 *)
NormalBoilingPoint::"usage"="NormalBoilingPoint[l] returns
the boiling point of liquid l at normal atmospheric pressure.";
NormalBoilingPoint[He]=4.222;
NormalBoilingPoint[H]=20.28;
NormalBoilingPoint[Ne]=27.09;
NormalBoilingPoint[N]=77.36;
NormalBoilingPoint[O]=90.19;
NormalBoilingPoint[Ar]=87.28;
NormalBoilingPoint[Kr]=119.77;
NormalBoilingPoint[Xe]=165.04;
NormalBoilingPoint[CH3]=111.63;
NormalBoilingPoint[NH3]=239.81;

(* Latent heat of vaporisation, Oxford Instruments p. 62 *)
LatentHeat::"usage"="LatentHeat[l] returns the latent
heat of vaporisation of liquid l. [J/kg]";
LatentHeat[He]=20.9 10^3;
LatentHeat[N]=198 10^3;


(* Simulation of complex systems *)
(* An attempt at object oriented programming with Mathematica *)

KrioNew::"usage"="KrioNew[o] creates a new system and stores it
into variable o. Other Krio function should be called with this
same variable o as the first argument.";
KrioNew[object_]:=Module[{},
object=Krio[{},{},{}];
];

SetAttributes[KrioNew, HoldFirst];

KrioGetSystemSize::"usage"="KrioGetSystemSize[o] returns the
number of parts in the system o.";
KrioGetSystemSize[object_]:=Length[ object[[1]] ];

KrioAddPart::"usage"="KrioAddPart[o, name, desc, temp:x] 
adds part to the system o. If temperature is defined, the part acts
like a heat reservoar with fixed temperature, i.e. it is either 
heat sink or hear source. Other parts will be in equilibrium,
i.e. the sum of all incoming and outgoing currents should be equal
to zero.";
KrioAddPart[object_, name_, description_, temperature_:x]:=Module[{parts,p},
parts=KrioGetSystemSize[object];
p=If[temperature===x, 
  KrioPart[name, description, T[parts+1]],
  KrioPart[name, description, temperature]
];
object[[1]]=Append[ReleaseHold[object][[1]], p];
] /; Head[ReleaseHold[object]]===Krio;

SetAttributes[KrioAddPart, HoldFirst];

KrioGetIndex::"usage"="KrioGetIndex[o, name] returns the index of
part name in the system.";
KrioGetIndex[object_, name_]:=Module[{index},
  index=Position[object[[1]] [[All,1]], name, {1}, Heads->False] [[1]];
  If[Length[index]==1, index[[1]], 0]
];
  
KrioAddTransfer::"usage"="KrioAddTransfer[o, name1, name2, desc, function]
adds a heat transfer process between two parts.
Function should take two temperatures as parameters and return
the heat current from part 1 to part 2.";
KrioAddTransfer[object_, name1_, name2_, desc_, function_Function]:=
Module[{index1, index2, temps},
  temps=object [[1]] [[All,3]];
  transforms=Table[ t[i] -> temps[[i]], {i,Length[temps]}];
  index1=KrioGetIndex[object, name1];
  index2=KrioGetIndex[object, name2];
  object[[2]]=Append[ReleaseHold[object][[2]], 
  KrioTransfer[J[index1, index2], 
	       function[index1, index2] /. transforms,
	       desc]
  ];
] /; Head[ReleaseHold[object]]===Krio && 
KrioGetIndex[object, name1]!=0 && KrioGetIndex[object, name2]!=0;

SetAttributes[KrioAddTransfer, HoldFirst];

KrioGetCurrents::"usage"="KrioGetCurrents[o] returns a list of all 
unknown currents.";
KrioGetCurrents[object_]:=Module[{},
  object [[2]] [[All,1]]
];

KrioGetTemperatures::"usage"="KrioGetTemperatures[o] returns a list of 
all unknown temperatures.";
KrioGetTemperatures[object_]:=Module[{},
  object [[1]] [[All,3]] ~ Cases ~ _T
];

KrioGetUnknowns::"usage"="KrioGetUnknowns[o] returns of all unknowns
in the given problem.";
KrioGetUnknowns[object_]:=Module[{},
  Union[KrioGetCurrents[object], KrioGetTemperatures[object]]
];

MAXCURRENT=100; (* Maximum currents for FindRoot *)
MINTEMP=4; (* Minimal temperature for FindRoot *)
MAXTEMP=300; (* Maximal temperature for FindRoot *)

KrioGetUnknownsFindRoot::"usage"="KrioGetUnknownsFindRoot[o]
returns a list of all unknowns in the given problem with
intervals of acceptable solutions. This form is suitable
for FindRoot function.";
KrioGetUnknownsFindRoot[object_]:=Module[
  {currents, temperatures, lenc, lent},
  currents=Union[KrioGetCurrents[object]];
  temperatures=Union[KrioGetTemperatures[object]];
  lenc=Length[currents];
  currents={currents,Table[0,{lenc}],Table[-MAXCURRENT,{lenc}],
  Table[+MAXCURRENT,{lenc}]}//Transpose;
  lent=Length[temperatures];
  temperatures={temperatures,Table[10,{lent}],Table[MINTEMP,{lent}],
  Table[MAXTEMP,{lent}]}//Transpose;
  Join[currents,temperatures]
];

KrioGetEquations::"usage"="KrioGetEquations[o] returns
the list of equations which need to be solved to fully
determine the thermal characteristics of the system.
We get both equations that describe the heat currents
as well as statements of the condition of the thermal
equilibrium.";
KrioGetEquations[object_]:=Module[
  {lhs, rhs, ndx, eqs1, eqs2, unknownT, in, out, 
  temps},

  (* Heat transfer *)
  lhs=Union[KrioGetCurrents[object]];
  Do[ ndx[ lhs[[i]] ]=i, {i,Length[lhs]}];
  rhs=Table[0, {Length[lhs]}];
  Map[ rhs[[ ndx[ #[[1]] ] ]] += #[[2]]&, object [[2]] ];
  eqs1=Equal[lhs, rhs] // Thread;

  (* Conditions of equilibrium *)
  unknownT=KrioGetTemperatures[object] [[All,1]];
  in=Map[ Plus@@Cases[lhs, J[_,#]]&, unknownT ];
  out=Map[ Plus@@Cases[lhs, J[#,_]]&, unknownT ];
  eqs2=Equal[in, out] // Thread;
  Join[eqs1, eqs2]
];

(* Helper functions for KrioAddTransfer *)

RadiationPlane::"usage"="RadiationPlane[eps1, eps2, A] returns
a function for KrioAddTransfer that describes thermal radiation
between plane surfaces of area A and emissivities eps1 and eps2.";
RadiationPlane[eps1_, eps2_, A_]:=Module[{},
HeatRadiationPlane[A, eps2, eps1, t[#2], t[#1]] &
];

RadiationCylinders::"usage"="RadiationCylinders[eps1, eps2, A1, A2]
returns a function for KrioAddTransfer that describes thermal
radiation between cylindrical surfaces of areas A1 and A2 with
emissivities eps1 and eps2.";
RadiationCylinders[eps1_, eps2_, A1_, A2_]:=Module[{},
HeatRadiationCylinders[A2, A1, eps2, eps1, t[#2], t[#1]] &
];

ConductionRod::"usage"="ConductionRod[m, l, r, n:1] returns
a function for KrioAddTransfer that describes heat conduction
through a circular rod made of material m of length l and
radius r. There can be n such rods.";
ConductionRod[material_, length_, radius_, number_:1]:=Module[{},
number ConductionHeatTransfer[GeomRod[radius, length],material,{t[#1],t[#2]}]&
];

ConductionPipe::"usage"="ConductionPipe[m, l, od, id, n:1] returns
a function for KrioAddTransfer that describes heat conduction
through a circular hollow pipe made of material m of length l, outer 
diameter od and inner diameter id. There can be n such pipes.";
ConductionPipe[material_, length_, od_, id_, number_:1]:=Module[{},
number ConductionHeatTransfer[GeomPipe[od, id, length],material,{t[#1],t[#2]}]&
];

KrioHeatTo[object_, to_]:=Module[{index,c},
  index=KrioGetIndex[object, to];
  c=KrioGetCurrents[object];
  (* TO DO *)
];

KrioSolve::"usage"="KrioSolve[o] solves the problem of thermal
equilibrium in system o. Function returns the replacement rules for
all the unknowns. Solution is also stored in the object o and
is used in KrioShowCurrents, KrioShowTemperatures, KrioCurrentAnalysis
and KrioHoldingTime.";
KrioSolve[object_]:=Module[{unk, sys, sol},
  unk=KrioGetUnknownsFindRoot[object];
  sys=KrioGetEquations[object];
  sol=FindRoot[sys, unk]//Chop;
  object [[3]]=sol;
  sol
];

SetAttributes[KrioSolve, HoldFirst];

KrioShowCurrents::"usage"="KrioShowCurrents[o] returns a matrix
with net thermal currents between parts of the system.";
KrioShowCurrents[object_]:=Module[
  {parts,np,allcurrents,unknowncurrents,othercurrents,currents},
  parts=object [[1]] [[All,2]];
  np=Length[parts];
  allcurrents=Table[J[i,j],{i,np},{j,np}]//Flatten;
  unknowncurrents=KrioGetCurrents[object];
  othercurrents=Complement[allcurrents, unknowncurrents];
  currents=Table[J[i,j],{i,np},{j,np}]/.Map[#->""&, othercurrents];
  (BlockMatrix[{ { {{""}}, {parts} }, { Transpose[{parts}], currents
  }}] /. object [[3]] )//MatrixForm
];

KrioShowTemperatures::"usage"="KrioShowTemperatures[o] returns
a table of temperatures of all parts of the system.";
KrioShowTemperatures[object_]:=Module[{},
  (Map[Apply[List,#]&,object [[1]]] /. object [[3]])//MatrixForm
];

KrioCurrentAnalysis::"usage"="KrioCurrentAnalysis[o, o] returns
a table with incoming and outcoming thermal currents for part p.
Also returnes is the net current to the part.";
KrioCurrentAnalysis[object_, partname_]:=Module[
{indeks, currents, currentsin, currentsout, totalin, totalout},
  indeks=KrioGetIndex[object, partname];
  currents=object [[2]];
  currentsin=Select[currents, #[[1,2]] == indeks &];
  currentsin=Map[ {#[[3]],#[[2]]}&, currentsin ];
  If[Length[currentsin]==0, currentsin={{0,0}}];
  totalin=Plus@@Column[currentsin,2];
  currentsout=Select[currents, #[[1,1]] == indeks &];
  currentsout=Map[ {#[[3]],#[[2]]}&, currentsout ];
  If[Length[currentsout]==0, currentsout={{0,0}}];
  totalout=Plus@@Column[currentsout,2];
  {"Currents in"->(currentsin//MatrixForm),
  "Total in"->totalin,
  "Currents out"->(currentsout//MatrixForm),
  "Total out"->totalout,
  "Net current"->(totalin-totalout)//Chop} /. object [[3]]
];

KrioHoldingTime::"usage"="KrioHoldingTime[o, p, l, c] returns
holding time for rezervoar p, filled with volume c of cryogenic liquid l.";
KrioHoldingTime[object_, partname_, liquid_, capacity_]:=Module[{},
  ca="Net current"/.KrioCurrentAnalysis[object, partname];
  capacity/(ca AmountBy1W[liquid]) /. object [[3]]
];

ShieldSolve::"usage"="ShieldSolve[l, d, t, l, eps, Tenv, Tsink] returns
a temperature profile on a cylindrical radiation shield of length l,
diameter d and thickness t. Thermal conductivity is l, emissivity
of the surface is eps, Tenv is the temperature of the environment,
while Tsink is the temperature where the shield is thermally sunk.";
ShieldSolve[length_, diameter_, thickness_, lambda_, eps_, Tenv_, Tsink_]:=
Module[{eq, shoot, target,sol, eq1, eq2},
  eq1=lambda thickness T''[z]
    ==eps StefanConstant (T[z]^4-Tenv^4); (* Pi diameter cancels out *)
  eq2=lambda Pi diameter thickness T'[0]
    ==Pi (diameter/2)^2 StefanConstant (t^4-Tenv^4);
  eq={eq1,eq2,T[0]==t};
  shoot[x_?NumberQ] := NDSolve[eq /. t -> x, T[z], {z, 0, length}][[1]];
  target[t_?NumberQ] := T[z] /. shoot[t] /. z -> length;
  (* Solve using the shooting method *)
  sol=FindRoot[target[t]==Tsink,{t, 0., Tenv},MaxIterations->1000];
  {sol,shoot[t]/.sol}
];

R0=0.0001; (* Trick to avoid problems at r=0 *)

RoundPlateSolve::"usage"="RoundPlateSolve[d, t, l, eps, Tenv, Tsink]
returns a temperature profile on a cylindrical plane of diameter d,
thickness t, thermal conducitvity l, emissivity eps. Temperature
Tenv is the environmental temperature and Tsink is the temperature
at the outside diameter of the plate where it is thermally sunk.";
RoundPlateSolve[diameter_, thickness_, lambda_, eps_, Tenv_, Tsink_]:=
Module[{eq, shoot, target, sol, eq1, eq2},
  eq1=lambda thickness (T'[r]+r T''[r])
    ==r eps StefanConstant (T[r]^4-Tenv^4); (* 2 Pi cancels out *)
  eq2=T'[R0]==0; 
  eq={eq1,eq2,T[R0]==t};
  shoot[x_?NumberQ] := NDSolve[eq /. t -> x, T[r], {r, R0, diameter/2}][[1]];
  target[t_?NumberQ] := T[r] /. shoot[t] /. r -> diameter/2;
  (* Solve using the shooting method *)
  sol=FindRoot[target[t]==Tsink,{t, 0., Tenv},MaxIterations->1000];
  {sol,shoot[t]/.sol}
];

EndPackage[];
