THESIS

J

 

lllllllllllllllllllIlllllllllllllllllll

3 1293 00892 7273

This is to certify that the

thesis entitled

AN INVERSE APPROACH TO THE ESTIMATION OF THE TISSUE
THERMAL PROPERTIES AND THE DETERMINATION OF THE OPTIMAL
TREATMENT TIME IN CRYOSURGICAL APPLICATIONS

presented by

Leslie Ann Scott

has been accepted towards fulfillment
of the requirements for

Master of Science degree in Mechanical Engineering

%mfleswr

0-7639 MS U is an Affirmative Action/Equal Opportunity Institution

 

Date July 30, 1993

 

 

 

LIBRARY
I Michigan State
UnIversIty

 

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

DATE DUE DATE DUE DATE DUE

 

 

 

‘x
:. i. ‘ -

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU Is An Affirmative Action/Equal Opportunity Institution
ammuna-n.‘

AN INVERSE APPROACH To THE ESTIMATION OF THE TISSUE
THERMAL PROPERTIES AND THE DETERMINATION OF THE OPTIMAL
TREATMENT TIME IN CRYOSURGICAL APPLICATIONS
By

Leslie Ann Scott

A THESIS
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Mechanical Engineering
1993

ABSTRACT
AN INVERSE APPROACH TO THE ESTIMATION OF THE TISSUE
THERMAL PROPERTIES AND THE DETERMINATION OF THE OPTIMAL
TREATMENT TIME IN CRYOSURGICAL APPLICATIONS
By

Leslie Ann Scott

Cryosurgery is a method Of destroying undesirable biological tissue, particularly cancerous
tumors, by freezing. Accurate estimation Of the thermal properties Of the tissue being frozen is
Often difficult due to its complex structure. Optimal duration Of treatment is Of extreme
importance, not only to ensure destruction Of the diseased tissue, but also to minimize loss Of the
surrounding healthy tissue. Methodologies are presented for the estimation Of the thermal
properties Of the tissue and the determination Of the Optimal treatment time. An infinite
homogeneous medium with constant thermal properties subject to a point heat sink was
considered. One-dimensional analytical solutions for dimensionless temperature in the frozen and
unfrozen regions were Obtained. Using these solutions, simulated data with added random errors
were used to evaluate me procedures. The estimated thermal properties using simulated data were
foundtobeinexeeflemagreememwimmepmperdesusedtogeneratethedata. Inthe
determination Of the Optimal treatment time, there was also excellent agreement between the

detenninedueannemtimeandmetimeusedtogeneratemesimulateddata

Copyright by
Leslie Ann Scott

1993

Dedicated with love tO
the memory Of my parents

Paul and Dorothy Scott

ACKNOWLEDGMENTS

I would like to express my eternal gratitude to the following people .....

TO my Godmother, Patricia Akin, for all Of her love, support and prayers.

TO my dear friend Mary Hawley, for being my ”big sister” and for always being there for me.
TO Connie Mosher, I couldn’t ask for a better friend.

To Julie Dickerson, for being such a good friend and for watching over me for Dorothy.

TO Doug Dotson, for always making me laugh (well almost always).

TO Barbara Penman. for all Of her pride and encouragement.

TO the Actons: Bob, Marge, Donna and Marilyn, for making me a part Of the family.

To Scott and Robin Hicks, for always being there to have fun.

TO Karen Calhoun, although She is no longer with us, she taught me the true meaning Of courage.
Her friendship and love will always remain in my bean.

TO Jim and Suzy locca. for being such great friends.
To Glenna Hackett, for all Of her understanding, friendship and support.

To Elaine Scott. for not only being my major professor and advisor, but for also being my friend.
i couldn’t have done this without her help and guidance.

To all Of my Engineering School friends: Laura. Doug, Debbie, Tammy, Roy, Rocky, Dave,
Amy, Kevin, Craig, Dan, and especially Kim Maier for being such a wonderful friend.

TO the members Of my advising committee: Professors James Beck, Elaine Scott. and Craig
Somerton.

TO The Dupont Company. for their financial support, and to the Department Of Mechanical
Engineering at MSU, both for their financial support and their encouragement.

To my life-long friend Jill Shannon, for the greatest gift I have ever received, my Goddaughter.

And most Of all to my Goddaughter, Casey Jill Shannon, for being the light Of my life.

TABLE OF CONTENTS

LIST OF TABLES

LIST OF FIGURES

NOMENCLATURE

CHAPTER 1
1.1

CHAPTER 2
2.1

2.2

2.3

2.4

2.5

CHAPTER3

3.1

INTRODUCTION
Objectives
LITERATURE REVIEW
Mechanisms Of Cell Destruction
Modeling the Heat Transfer in Living Tissues
2.2.1 Analytical Solutions to the Describing Differential Equations

2.2.2 Numerical Methods for Solving the Describing Differential
Equations

Experimental Investigations
Monitoring the Freezing Pro'cess
2.4.1 Ultrasound Imaging in Cryosurgery

2.4.2 An Alternative Method for Determining the Interface
Location

Optimization Of the Freezing Process
THEORETICAL METHODS

The Parameter Estimation and Inverse Problem Solution Techniques

vi

xii

10

10

11
12

13

3.2

3.3

3.4

CHAPTER4

4.1

4.2

4.3

4.4

vii
3.1.1 The Box-Kanemasu Interpolation Method
Mathematical Description of the Direct Problem
Estimation Of the Thermal Properties

3.3.1 Determination of the Sensitivity Coefficients

The Inverse Problem of Determining the Optimal Cryosurgical
Treatment Time

ANALYTICAL PROCEDURE
The Sensitivity Coefficient Analysis for the Thermal Properties
4.1.1 Magnitudes of the Sensitivity Coeflicients
4.1.2 Linear Dependence Between Sensitivity Coefficients
The SumoquuaresPunctionforthe'I‘bermalProperties
Estimation of the Thermal Properties
4.3.1 Input for the Box-Kanemanr Interpolation Method
4.3.2 Individual Estimation of the Thermal Properties
4.3.3 Simultaneous Estimation of the Thermal Properties

‘I‘heDeterminationProcedurefortheOptimal Cryosurgical
TreatmentTime

4.4.1 Sensitivity Coefficient Analysis for the Optimal
Cryosurgical Treatment Time

4.4.2 The Sum of Squares Function for the Optimal
Cryosurgical TreatmentTime

4.4.3 Modification of the Box-Kanemasu Interpolation
Method Program

4.4.4 Input for the Box-Kanemasu Interpolation Method
4.4.5 Determination Of the Optimal Cryosurgical Treatment Time
4.4.6 Use of Prior Information Obtained from a Different Radius

4.4.7 Use of Prior Information Obtained from two Different Radii

16
18
23

23

27
31
31
33

37

39

39

41

42

42

45

3

3

47

49

49

viii

CHAPTER 5 RESULTS AND DISCUSSION 51
5.1 Estimation of the Thermal Properties 51
5.1.1 Individual Estimation of the Thermal Properties 51
5.1.2 Simultaneous Estimation of the Thermal Properties 57
5.1.3 Comparison Of Individual and Simultaneous Estimates 63
5.2 The Optimal Cryosurgical Treatment Time 63
5.2.1 Determination of the Optimal Cryosurgical Treatment Time
using Prior Information from the Same Radius 63
5.2.2 Determination of the Optimal Cryosurgical Treatment Time

using Prior Information from a Different Radius 68
CHAPTER 6 SUMMARY AND CONCLUSIONS 73
6.1 Estimation Of the Thermal Properties 73
6.2 Determination Of the Optimal Cryosurgical Treatment Time 74
CHAPTER 7 RECOMMENDATIONS FOR FUTURE WORK 76
APPENDIX A TI-E FORTRAN PROGRAM SEN SE.FOR 78
APPENDIX B TI-E FORTRAN PROGRAM SQUARESFOR 86
APPENDIX C TI-E FORTRAN PROGRAM NLINA.FOR 91
APPENDIX D TIE FORTRAN PROGRAM MOD.FOR 112

APPENDIX E THE SUBROUTINES MODEL AND SENSE FROM NLINA.FOR

MODIFED FOR THE DETERMINATION OF TIE OPTIMAL

TREATMENT TIME 125
APPENDIX F TI-E FORTRAN PROGRAM MODC.FOR 134
APPENDIX G TIE FORTRAN PROGRAM RAD.FOR 140

APPENDIX B TIE SUBROU'TINES MODEL AND SENSE FROM NLINA.FOR

BIBLIOGRAPHY

MODIFIED FOR THE DETERMINATION OF THE OPTIMAL
TREATMENT TIME WITH PRIOR INFORMATION FROM
A DIFFERENT RADIUS 145

161

Table 5.1.

Table 5.2.

Table 5.3.

Table 5.4.

Table 5.5.

Table 5.6.

Table 5.7.

Table 5.8.

Table 5.9.

Table 5.10.

Table 5.11.

Table 5.12.

Table 5.13.

Table 5.14.

Table 5.15.

LIST OF TABLES

EsfimafionofdreDimensionlessLaeruHeatofFusimL‘

Estimation of the Dimensionless Thermal Conductivity, k,,

Estimation of the Dimensionless Thermal Diffusivity, (1,,

Number of Iterations Required for Convergence in the Estimation Of L‘
Number Of Iterations Required for Convergence in the Estimation Of lg,
Number of Iterations Required for Convergence in the Estimation of (1,,

Simultaneous Estimation Of the Dimensionless Latent Heat Of Fusion, L‘,
and the Dimensionless Thermal Diffusivity. (1,,

Simultaneous Estimation Of the Dimensionless Thermal Conductivity. k...
and the Dimemionless Thermal Diffusivity, a,

Number of Iterations Required for Convergence in the Simultaneous
Estimation of L’ and a.

Number of Iterations Required for Convergence in the Simultaneous
Estimation of k, and (1,,

Estimation Ofthe Optimal Treatment Time, rausing SimulatedData
frmnboththeFrozenmdUnfrozenRegions

Estimation of the Optimal Treatment Timc. 1,, using Simulated Data
from the Frozen Region Only

Estimation of the Optimal Treatment Time, r,, using Simulated Data
from the Unfrozen Region Only

Number of Iteratitms Required for Convergence in the Estimation of re
usingSimulatedDatafromboththeFrozenandUnfrozenRegions

Number Of Iterations Required for Convergence in the Estimation of re
using Simulated Data from the Frozen Region Only

ix

52

52

53

56

56

57

58

59

62

62

65

67

68

Table 5.16.

Table 5.17.

Table 5.18.

Table 5.19.

Table 5.20.

I

Number of Iterations Required for Convergence in the Estimation of I,
using Simulated Data from the Unfrozen Region Only

Estimation Of the Optimal Treatment Time. 1‘, using Prior Information
Obtained from a Different Radius

Estimation of the Optimal Treatment Time. I, using Prior Information
Obtained from Two Diflerent Radii

Number Of Iterations Required for Convergence in the Estimation of I,
using Prior Information Obtained from a Different Radius

Number Of Iterations Required for Convergence in the Estimation of I,
using Prior Information Obtained from Two Different Radii

68

69

7O

71

71

Figure 3.1.

Figure 4.1.
Figure 4.2.
Figure 4.3.
Figure 4.4.
Figure 4.5.

Figure 4.6.

Figure 4.7.

Figure 4.8.

Figure 4.9.

LIST OF FIGURES

Propagation Of the Freezing Front, :0). as a Function ofTime
and Radius in an Infinite Homogeneous Medium with
Constam Thermal Properties

Dimensionless Sensitivity Coefficients versus 11
Dimensionless Sensitivity Coefficients X 1:. versus X L.
Dimensionless Sensitivity Coefficients X9 versus X L.
Dimensionless Sensitivity Coefficients Xe, versus X“
‘I‘heSumoquuaresFunctionversusL'

T'heSumoquuares Functionvernrskdandau
Dimensionless Temperature of the Cryosurgically Frozen Tissue
versus Time at a Radius Of 0.1 Meters, Given a Treatment
Time I, a 0.185 Seconds

Dimensionless Sensitivity Coefficient X5 versus 11

The Sum of Squares Function versus rc

14

33

35

35

36
38

38

t

45

NOMENCLATURE

Coefficient used in the assumed solution Of the differential equation describing the
temperature of the frozen region

Scalar used in the Box-Kanemasu Interpolation Method

Coefficient used in the assrmed solution Of the differential equation describing the
temperature Of the fiozen region

Vector of estimated parameter values

Coefficient used in the assumed solution of the differential equation describing the
temperature Of the unfiozen region

Specific heat of blood (kJ/kg-K)
Specific heat of the frozen tissue (kJ/kgK)

Coefficient used in the assumed solution of the differential equation describing the
temperature of the unfrozen region

Exponential fimction

Complimentary error function

Transcendental equation for the dimensionless freezing fiont location A.
Scalar used in the Box-Mam Interpolation Method

Scalar interpolation factor

Thermal conductivity of the Imfiozen tissue (W/m-K)

Thermal conductivity of the frozen tisnre (W/mK)

Dimensionless thermal conductivity

IaterrtheatoffusionofthetissueGJ/kg)

q-s

Dimensionless latent heat Of fusion

I-Ieat sink (W/m)

Dimensionless heat sink coefficient

Metabolic heat generation rate (kW/m3)

Radius (m)

Sum Of squares function

Freezing from location (m)

Vector of calculated data

Systemic arterial blood temperature (°C)

Initial temperature of the tissue prior to freezing (°C)
Temperature in the unfrozen region (‘0
Temperature at the interface between the fiozen and unfrozen regions (°C)
Temperature in the frozen region (°C)

Time (s)

Optimal uearment cooling time (s)

Actual time minimum temperature is achieved due to difiusion (s)
Blood perfusion rate (kg/m’s)

Sensitivity coeflicient matrix

Dimensionless sensitivity coemcient

Dimensionless sensitivity coefficient for lg,
Dimensionless sensitivity coefficient for L'
Dimensionless sensitivity coemcient for rc
Dimensionless semitivity coefficient for 0,,

Vector of measurement data

xiv
Process variable

Infinity

Scalar used in the Box-Kanemasu Interpolation Method
Thermal diffusivity of the unfiozen tissue (mzls)
Thermal diffusivity of the frozen tissue (m’ls)
Dimensionless thermal diffusivity

Parameter to be estimated

Parameter vector

Vector used in the Box-Kanemasu Interpolation Method
Dimensionless similarity variable

Dimensionless similarity variable

Dimensionless temperature in the unfrozen region
Dimensionless temperature in the frozen region
Dimensionless freezing from location

All parameters not being estimated

Density Of the tissue (kg/m3)

Standard deviation

Arterial
Blood
Cryosurgical
Initial

Unfrozen (liquid)

Melting

mg Metabolic generation

min Minimum

p Pressure (constant)

3 Frozen (solid)

s1 Indicates dimensionless ratio of solid to liquid
W

k Iteration number

T TmsPose

* Indicates dimensionless quantity

CHAPTER 1

INTRODUCTION

Cryosurgery is a medical technique Of destroying undesirable biological tissue, particularly
cancerous tumors, by freezing tO very low temperatures (typically around -50°C). The use Of
freezing in the treatment of malignancy dates back to the 1850s when iced saline solutions were
used to treat advanced cases Of breast and uterine cervical cancer (Gage, 1992). It was found that
the freezing treatment provided a relief Of pain, a reduction in tumor size, and a decrease in
bleeding and discharge. The availability Of liquified gases in the late 18005 further developed the
field Of freezing therapy. Liquid air was often used to treat a variety Of skin disorders, including
skin cancer, with favorable results reported.

Pioneering work in this area continued during the years 1936 to 1940 (Gage, 1992).
Patients with large incurable cancers of the breast and uterine cervix were treated with irrigations
of iced saline solution through hollow instruments ill contact with the tumor. The benefits again
included a relief Of pain and a reduction Of tumor size, just as reported nearly a century earlier.
Medicalresearchinthisaleawas interruptedby World WarIIanddiantcontinue forseveral
years after due to the association in scientists’ minds with the Nazi’s infamous hypothermia
experiments in concentration camps (Rubinsky, 1986b).

Research resumed in the 1950s with impressive results reported, but it wasn’t until 1961
that modern cryosurgery was born. This was due to the development Of a new cryosurgical

apparatus by Irving Cooper and AS. Lee. This device, the cryoprobe, is a small hollow cylinder

2
that is vacuum insulated everywhere except the tip. A cryogen, usually liquid nitrogen. is
circulated through the cryoprobe. Freezing of the undesirable tissue is achieved through the
placement Of the cryoprobe either in contact with the surface Of the tumor or directly into the
mmorbyplmcmm,resulfingindreremovalofheatfiommesunoundingfissue. Asthisheatis
removed, afieezing front propagates outward fiom theprobe resulting indestruction Ofthe tissue.

Cryosurgery Offers many advantages over traditional forms Of treatment Of cancer. It
does not require the resection of large volumes Of sunounding healthy tissue. Often, anesthesia
and surgery are not necessary. Also, because fieezing can be localized, multiple lesions within
an organ can be heated individually. However, for a cryosurgical procedure to be successful, it
isextremely importanttobeabletopredicttheareaoftissuenecrosis. TOdO so requiresthatthe
thermalconditionsthatpromotethedestructionofthetissue,aswellasthethumalhistoryand
extent Of fieezing during the procedure, are well understood (Rubinsky, 1986b). It is essential that
destruction of the entire tumor is achieved while damage to the sunounding healthy tissue is
minimized. Therefore, the propagation Of the freezing fiont must be precisely determined during
the freezing process.

Themostprominentuseofcryosmgeryisinthetreannentofskincancer. Inthiscase,
the location Of the freezing front is easily seen by the surgeon. Other applications include the
treatment Of prostatic and uterine cervical cancers, the ueatment Of Parkinson’s disease by
destroyinglesionswitbintbebrain.andthedestructionoftumorswithintheliver. Thesedeep
bodylocafionsdonotpermhmevisualizafionofurefieezingfiombymesurgeon The
heterogeneous nature and large blood supply (typical of malignant tissue) of the tumor further
complicatethefreezingprocess. Inmanycases,thedurationoftreatmentandamormtoftissue
destroyed has become almost an art form for the surgeon. with experience gained through
previously performed procedures.

TO accurately determine the optimal duration of cryosurgical treatment, specific for each

3
size and type Of tissue to be destroyed, appropriate mathematical models are needed to describe

the heat transfer within the tissue. The use of mathematical models also requires the accurate
knowledge Of the thermal properties of the tissue to be frozen. however, the values Of these tissue
thermal properties are Often unavailable. Therefore, the development Of methods that allow for
fluaccurateesfimafionofdwfissuednrmflpmperfiesandopfimalmrmfimofueamemis
essential.

In this investigation, mathematical models were chosen to describe the heat transfer within
the tissue undergoing cryosurgical freezing. Simplifying assumptions were made to allow for the
attainment Of analytical solutions to these describing differential equations. A minimization
procedure, the Box-Kanemasu Interpolation Method, was used to estimate the thermal properties
and to determine the Optimal duration Of treatment. Although the problem was simplified. the
primary goal Of this study was to test the methodologies for accuracy and reliance.

1.1 Objectives

The Objective of this investigation was two fold. First, the minimization prowdure was
used to estimate the tissue thermal properties. Specifically, the properties to be estimated were
the latent heat Of fusion, thermal conductivity, and thermal diffusivity. Second, the minimization
pmcedumwasusedmdeterminedreopfimalueaunemfimemquimdmachieveadesimd
minimum temperature at a specified radius location of the tumor. In both cases, simulated
measurement data were required as input-for the procedure. This allowed for the comparison
betweentheestimatedvaluesandtheactualvaluesusedtogeneratethesimulamddata.

In the sections that follow, a literature review is presented in Chapter 2. In Orapter 3.
Theoretical Methods, the Box-Kanemasu Interpolation Method is described. with all necessary
mathematical expressions introduced. The methods used to estimate the thennal properties and
todeterminetheoptimaltreatmenttimearepresentedinalapter4, AnalyticalProcedureS. The

results ofthisinvestigationarepresentedanddiscussedinalapter5,withtheconclusionsgiven

4
in Chapter 6, Summary and Conclusiom. Flnally, Chapter 7 provides recommendations for future

work in this area. All computer programs required for this study are located in the Appendices.

CHAPTER 2

LITERATURE REVIEW

Cryosurgery is a medical treatment in which malignant and other biological tissue is
destroyed by freezing to very low temperatures. It is essential that the diseased tissue is destroyed
while maintaining as much healthy tissue as possible. The outcome Of the cryosurgical treatment

is dependent upon the cooling and warming rates imposed at the freezing/thawing fiont.

2.1 Mechanisms of Cell Destruction

The cells that comprise the undesirable tissue are destroyed by two distinctly different
methods, depending upon the rate of cooling/warming. When cells are subject to a slow rate of
cooling (10°C/min). they often remain unfrozen. yet supercooled (Savic, 1984). Ice forms in the
vasculature system first while the cells adjacent to the blood vessels remain unfrozen (Rubinsky
and Eto, 1989). Since the supercooled water has a higher vapor pressure than the ice crystals, the
cell equilibrates by losing water, resulting in dehydration Of the cell. There is a subsequent
concentration ofthe solutes within the cell that. when high enough, leads to the death ofthe cell.
Inaddifiomwbentherateofcoolingissiow,theformationoficeinthevasculaulre realltsin
expansionofthevessels. Thiscausesalossofstrucurraiintegrityofthevessel,whichmaynot
be functional when thawed.

When rapid cooling occurs (100°Clmin), the cell is unable to equilibrate quickly enough

6
through dehydration. Therefore, equilibrium is achieved by the formation of intracellular ice. The

small ice crystals produced through rapid cooling are likely to fuse with other small crystals
withinthecell. Thisresultsindrenrphueofthecehmembraneandthedeathofthecell. While
a slow rate of cooling may not always result in the death of cells, as in fiostbite injuries.
intracellular ice formation is almost invariably lethal to cells (Comini and Del Guidice, 1976).
Smaflaystalspmducedbympidcoolhghaveatendencymrecrystalfizedunnguc
thawing process, especially if the warming is slow. This results in the breaking of more cell
membranes, killing more cells. Therefore, more cells are killed during a slow thawing process.
In cryosurgical applications, the ideal case is rapid freezing and slow drawing to achieve the

maximum percentage of cell death.

2.2 Modeling the Heat Transfer in Living Tissues

The complicated and only partially understood mechanisms involved in cell freezing and
destnlction make the modeling of the heat transfer within the tisme quite difficult. Another
obstacle is the nonhomogeneous nature of the tissue and the temperature dependent, relatively
unknown. thermal properties. Many attempts have been made to predict the temperature profiles
within the tissue undergoing fieezing.

2.2.1 Analytical Solutions to the Describing Differential Equations

Oneofthetraditional waystodeterminethetemperatureprofilesduringcryosurgeryis
to solve the treat conduction problem for a two-phase system with a moving boundary between
the frozen and unfrozen regiom. A major obstacle is the lack of accurate published data regarding
the tissue thermal properties. Also, factors petalliar to biological systems must be incorporated
huodedescfibingequaflomnwhasmemoodpefiusionandmembohcheagaerafimmes

(Filippov and Vasil’kov, 1979).

7

In work presented by Rubinsky and Shitzer (1976) the bio-heat eqration was used to
describe the heat transfer in the unfrozen region of the tissue. In the frozen region the diffusion
of heat was described by the heat equation. The assumptions of homogeneity and constant
thermal properties were made. Another assumption was that the blood perfusion and volumetric
metabolic rates were considered to remain constant throughout the cooling process prior to
freezing. Upon freezing, these quantities go to zero. The bio-heat equation describing the heat

transfer in the unfrozen region is

31 air MC. 4
= __ T - ——n
.37 0—2 + PC,( . T) + pCP (2.1)

where

w, = blood perfusion rate

C, = specific heat of blood

7', = systemic arterial blood temperanue

q.‘ = metabolic heat generation rate

Simplifying assumptions were made that allowed for the attainment of analytical solutiom
to both the bio-heat and treat transfer equations. The results of this investigation showed that the
temperature gradient at the freezing fiont becomes larger as the volumetric metabolic and the
blood perfusion rates are increased. thereby decreasing the velocity of the freezing front
propagation. Theeffectofthebioodperfusionratewas found tobemuch strongerthanthe
volumetric metabolic rate. It was concluded that q... could be neglected without causing
significant errors in the results. However. the blood perfusion effects were significant and could
not be ignored.

To predict the maximum size of the cryolesion (fiozen portion of the tissue), a steady state
version of the bio-heat transfer equation was presented by Cooper and Trezek (1972) to describe

thetemperaturefieldsirlbothmefiozenandunfiozenregionsofthetissm.subjecttodre

8
appropriateboundary conditions foreachregion. Analytical solutionswereobtainedandusedm
predict the steady state or maximum cryolesion size for different shapes of cryoprobes: planar.
cylindrical. and spherical. From previous work (Cooper and Trezek. 1970) they had also formd
thattheeffectsoftheheatcapacityofthetismewerenegligibleinboththefrozenandunfrozen
regions. However.thelatemheatoffusioneffectsresrufingfiomflrephasechmgewere
significant.

A semi-infinite slab. initially at a uniform temperature, was the considered geometry in
a freezing simulation performed by Warren et al. (1974). This flat geometry was formed by an
elastic bladderpressed againstthetissue. Thebladderwasthencooled intemallybyacryogen.
eitherliquid nitrogenorfreon. Integral equationswereusedtodescribetheheattransferinthe
frozen and unfrozen regions. The temperature profiles were approximated by polynomials subject
to either a temperature or convective boundary condition. This allowed for the attainment of
simple closed form expressions for both the tissue temperature and the frozen/unfrozen interface
position. They found that this technique was particularly valuable for demonstrating the
effectiveness of coolants.

Hrycaltet al. (1975) used the cylindrical form ofthe heattransferequationto describethe
freezing process during cryosurgery. Analytical solutions of the temperature distribution were
obtained and used in the Neumann’s solution to the frost penetration problem in a semi-infinite
slab. Fromthissolutimthetimerequiredtoachievefieezingatagivendepthwasobtainedand
couldbeusedasanestimateofthecryosurgicaltreatmenttime.

2.2.2 Numerical Methods for Solving the Describing Differential Equations

In an investigation performed by Hayes and Diller (1982), the bio-heat transfer equation
wasusedmdesaibemeheatwnducnonprocessmdnhummbodysubjeaedmexheme cold.
It was felt that this type of analysis could be used as a predictive tool in many applications.

including cryosurgery. A finite element model for a composite human was used to solve for the

9
temperature profile within the body. To avoid the complexity of a three-dimensional model of

a human. the model was simplified by using a two-dimensional. axially symmetric model. The
resultsofthis study showedthatmeeffectsoflatanheatreleaseandbloodperfiisionkeepthe
tissuewannerthaniftheywerenegleeted. Also.asthetissueisundergoingaphasechange,the
effects of latent heat greatly influence the predicted temperature field. It was concluded that the
latent heat had much more effect on the temperature field than the blood perfusion rate.

The finite element method was also uwd by Comini and Del Guidice (1976) to solve a
two-dimensional bio-heat transfer equation coupled with distributed convection at the exposed
surface of the tissue. Temperature profiles were determined for various applications. such as the
treatrnentoflarge angiomasininfants. Inthiscasethethermalpropertiesofthebrainandskull
had to be estimated. thereby re-affirming the need for accurate determination of the
thermophysical properties.

In work done by Budman et al. (1990), the heat transfer equation was used to describe
the conductive heat transfer in the frozen and unfrozen regions. Also. they assumed the presence
of an intermediate range of frozen plus unfrozen phases. resulting in a modified heat equation to
include the effects of latent heat release. Solutions to the three describing differential equations

were obtained through the use of the Runge-Kutta Method.

2.3 Experimental Investigations

In experimental work done by Cooper and Petrovic (1974) a solution of 1.5% gelatin.
98.5% waterwasusedasatestmediumtosimulatetissue. Alsoincludedinthismedium wasa
liquid crystal sheet. Thiscrystal sheethadthefeatmeofdisplayingbrilliantchanges incolorover
discrete temperature bands. Using liquid nitrogen as the cryogen and a cryoprobe to provide
fieezing. photographs were taken of the frozen region and the various isotherms displayed on the

10

crystalsheetatdifferenttimeintervals UsingananalyticalmodelproposedbyCooperand
Trezek that neglects heat capacity in both the frozen and unfrozen regions, Cooper and Petrovic
found that the rate of growth of the frozen region. determined experimentally. compared within
19%ofthosepredictedusingthetheoreticalmodel.

Thepmcessoffieezinginflssueisaffeaedbyflreuansponofwaterfiommeceusinto
the surrounding vasculature. Therefore, this mass transport of water must be considered in the
formulation of the energy equation. This mass transport process was experimentally studied by
Rubinsky and Etc (1989) with the use of a "Krogh Cylinder". This tmit consists of a cylindrical
blood vessel surrounded by tissue. The tissue is represented by a solution of electrolytes in water.

From this work an expression for the change in the radius of the blood vessel was determined.

2.4 Monitoring the Freezing Process

From experiments performed by Augustynowicz and Gage (1985) it was determined that
thetemperamreatadepthwithinthetissuewas always lowerthanthatatanequidistantsiteon
the surface. In general. clinicians have based their judgement on the depth of freezing to be
approximatelyequaltothelateralspreadoffiostfromtheprobe. Thisemphasizestheneedfor
wcumtemefingofflefieezingprocessmracwmmmemodsofdetemmingmeopfimal
duration of treatment for a given tumor size, in cryosurgical procedures.

2.4.1 Ultrasound Imaging in Cryosurgery

ItisabsolumlycmcialmproduceapredicmMemeaofmcmsismdwcryosurgical
treatment of cancer. Insufficient freezing leaves viable cancer cells while over freezing can have
disastrous consequences (Gilbert et al., 1984). Therefore. the growth of the freezing front must
be accurately determined during freezing. Previously. thermocouples inserted near the margins

ofthemmorhavebeenusedtomonitorthefreezingprocess. Thisallows onlyalimiwdnumber

l 1
of points to be monitored and is difficult when the tumor is heterogeneous or large blood vessels

are in the vicinity. Also. needles containing thermocouples cannot be easily imerted into tumors
locateddeepinthebody.suchasthebrainorliver.

Ultrasound, an imaging method which uses sound waves. has recently been used to
continuously monitor the position of the phase change interface during cryosurgery. Ultrasound
devices produce images by analyzing the acoustic energy reflected from tissue through which an
acoustic pulse has been transmitted. Gilbert et al. (1984) performed experiments on simulated
organ tissue (transparent bovine gelatin). They found that the frozen region produced during
cryosurgery is clearly visible under ultrasonic imaging. Therefore, while thermocouples only
permit a limited number of points to be monitored, which is inadequate for accurately predicting
the freezing process. ultrasound allows for continuous monitoring of the position of the phase
change interface.

Experiments performed on laboratory animals have shown that the frozen/unfrozen
interface is a strong reflector of acoustic energy (Rubinsky, 1986a). Experiments performed on
gelatin samples in both planar and hemispherical freezing processes have also shown that the
change of phase interface is easily identifiable by ultrasonic monitoring (Gilbert et al., 1985). It
has been determined that ultrasound can be used to continuously monitor the transiait position
of the frozen/unfrozen interface during cryosurgery.

2.4.2 An Alternative Method for Determining the Interface Location

Visualization techniques. particularly ultrasonic monitoring, allow relatively accurate
control of freezing. However, this monitoring technique provides only two-dimensional
information of a three-dirnennonal process. As an alternative. a microprocessor data collection
system was presented by Savic (1984). Since the exact mechanism of cryosurgery is very
complicated and not fully understood. it was decided to experimentally collect data on the

changing physical environment during freezing. This data was then compared with the results of

12
the biopsy of the tissue to find a relationship between the percentage of cells killed and the

parameters describing the physical environment of the tissue. thereby determining which
combinations of parameters accurately predict cell death. The phenomena being simultaneously
momwmdwemmefismetanpemmwusmgmemomuplesmdflmfismemsistmusmgamedle
probe. which measured the increase in resistance with a decrease in temperature. While this
method shows promise. considerable work is needed in this area
2.5 Optimization of the fleecing Process

Keanini andRubinsky (1992)presentedatechniquethatminimizesumrecessaryfreezing
byoptimizingtlrenumberandsizeofcryoprobesusedintheprocedure. Thistechniqueusedthe
Simplex algorithm, a minimization technique used to determine function minima. In this case the
functiontobeminimized wasthevolumeofhealthytissuedestroyedduringthefieezingprocess.
This function was assumed to depend on three independent variables: the number of probes, the
probe diameter. and the probe active length. Although this method shows promise. considerable
work is yettobedone. Forinstance,definingthe functiontobeminimizedisproblem specific.
depending upon the type of tissue to be frozen. Therefore. accurate biophysical and bio-heat

transfer models are needed, as are efficient minimization algorithms.

CHAPTER 3

THEORETICAL METHODS

There are two major aspects in this study of the cryosurgical freezing of undesirable
tissue. The first is the estimation of the thermal pr'Operties of the tissue. namely the latent heat
of fusion. thermal conductivity, and thermal diffusivity. From the Literature Review presented
in Chapter 2 it was found that the effects of the latent heat of fusion were quite significant in the
freezing process. The second aspect involves the determination of the optimal cryosurgical
treatment time required to achieve a desired minimum temperature at a specified location.

The actual problem of determining the thermal properties and optimal treatment time has
been simplified to test the estimation procedure for accuracy and reliance. In both analyses. the
tissue to be frozen was heated as an infinite homogeneous medium. initially at a uniform
temperature. Also. the effects of blood perfusion and metabolic heat generation were ignored.
Due to the assumption of a homogeneous medium. the cryosurgical probe was considered to be
located atthegeometric centerofthetissue and wasmodeled as apointheatsinkthatfreezesthe
surrormdingtissue. Thisresultsintwophaseswithinthemediumzafrozenregionandan
unfiozenregionsepmatedbyminterface,whichisknownasflrefieezingfiom. Thethermal
propertiesofthefiozenandtmfrozenregiom were assumedtobeconstantwithineachphasebut
different between phases. Therefore. the freezing front. so). propagates spherically outward from

the probe in a one-dimensional fashion as a fimction of time and radius as shown in Frgure 3.1.

13

l4

Banner—1395.28.23
3.3,. 250895: 82...... .a s 8.3.8. 25% SE. a .8» flan—gnu... 8‘98: .2 as»:

R

 

    

 

:23. 83:

 

8:3. 88:5

 

gm «3. 23..

\

 

 

823.... .8»: 528
5:6: «.388: 32:...

 

15
The determination of the temperature profiles of the frozen and unfrozen regions from the

initial. boundary. and interface conditions is considered to be mathematically well-posed and is
called a direct problem. The estimation of the thermal properties from internal temperature
measurement data is commonly referred to as parameter estimation. while the determination of
the optimal treatment time from internal temperature measurement data is called an inverse
problem. These problems are considered to be mathematically ill-posed. The one-dimensional
temperature profiles for the frozen and unfrozen regions. obtained from the direct problem. are
required for the estimation of the thermal properties and the determination of the optimal treatment
time. The assumptions of constant thermal properties and homogeneity of the tissue allowed for

the attainment of exact solutions to the describing differential equations.
3.1 The Parameter Estimation and Inverse Problem Solution Techniques

The estimation of the thermal properties of the tissue and the determination of the optimal
treatment time involved the use of the same minimization procedure in both analyses. This
involves minimizing an objective function, such as the sum of squares function. with respect to

a particular parameter of interest. The sum of squares function is

s = [Y - T(B)l’ [Y - 7(3)] (11)
where Yis avectoroftemperature measurementdata. and Tisavectorofcalculated temperature
values from the describing model and is a function ofthe unknown parameters contained in the
parameter vector B. In the parameter estimation problem, B contains the unknown thermal
properties.andintheinverseproblem. Bcontainstheoptimal treatrnenttime. Tominimizethe
sum of squares function. equation (3.1) is differentiated with respect to the unknown parameters

andsetequaltozero. Theresultingexpressicnis

16

V55 = ZI-X'(B)][Y ' 7(3)] = 0 (33)
The vector B contains the true parameter values. and X03) is the sensitivity coefficient matrix.

defined as
X(B) - [tartan]T (3.3)

Many times there exists prior information of the parameter to be estimated. To include

this prior information. the sum of squares ftmction is modified as follows:

s = [Y - raw [Y - 7(3)] + [B - 01' [B -b1 (3")
where B is a vector of the actual values of the parameters being estimated. obtained from prior
information. and b is a vector of the estimated values of the parameters.

3.1.1 The Box-Kanemasu Interpolation Method

Aprobls- afisesintheleastsquaresmethodwhenflremathemaficalmodelisnonlinear
in terms of the parameters. In this case. it is not possible to explicitly solve equation (3.2) for the
parameter vector B. The Gauss Method of Minimization (Beck and Arnold. 1977) is a simple and
effective method that provides a linear approximation to the nonlinear model. To transform
equation (3.2) into an iterative form. two approximations are used: 1) the sensitivity coefficient
matrix. X05). is replaced with X(b). where b is an estimate of B and 2) the vector of calculated

valueS. 7(5). is approximated by using a Taylor expansion of Tm) about 0 as follows:
1(a) - m) + [warm]T (B - b) + . . . . (35)
Neglecting the higher order terms of the Taylor series results in the following expression for v.5:
VnS 5 X'(b)[Y - 1'0) - X(b)(B - b)] E 0 (3.6)

whichinnowlinearintermsoftheparametervectorfl.
Because the Gauss Method uses a linear approximation of Ni). oscillations and

nonconvergence can sometimes occur in the iterative process. The Box-Kanemasu Interpolation

17
Method. also presented by Beck and Arnold (1977). is a modification of the Gauss Method that

eliminates this problem of nonconvergence. The Box-Kanemasu Method provides an iterative

procedure for the estimation of B with b as follows:

W) = rm + a “may” <17)
where

A,” = [X’(b)X(b)]"[X""(b)(Y - 1%)» (33)

The iteration number is k. and It“) is a scalar interpolation factor.

In the Box-Kanemasu Method. the sum of squares function is approximated at each

iteration using

S = a0 + alh + azh2 (3.9)
Theconstants ao.a,.anda2arecharacteristicofeachiterationandareequalto

a0 = 50‘". a1 = -20 ('0, a2 = is!” - s3" + 26 mania-2 (3-10 “-99
where

G“) = [A,W’]T[Xr(b)x(b)]wA,bm (3.11)

Initiallya=1.andSoandS¢_,arethevaluesofSwithh=Oandh=1respectively. This

approximate form of s is then minimized to calculate am
W” = G“’a’lS.“’ - 8.5" + 26%)" (3.12)

The calculated ha“) is then used in equation (3.7) for the (k+1)st iteration of b.

AcheckismadeaftereachiterationtoeonfirmthatSisindeeddecreasingbyensuring

S a?) < S on) (3.13)

with a being made sufficiently small for this to occur. Iteration proceeds until there is little

18

difference between D“) and b“).

In this investigation. the use of the Box-Kanemasu Interpolation Method required solutions
to the mathematical models describing the temperature profiles within the frozen and unfiozen
regions of the tissue being frozen. Also required were the sensitivity coefficients with respect to
thethermalpropertiestobeestimatedandtheoptimalcryosurgicaltreatrnenttimetobe

determined. as well as simulated experimental measurement data to be used as input.

3.2 Mathematical Description of The Direct Problem

The mathematical model used in the Box-Kanemasu Method is the same for both the
estimation of the thermal properties and the determination of the optimal treatment time. Since
the freezing process results in the presence of two phases within the medium, the problem must
be described with two equatiom: one for the frozen (solid) region and one for the unfrozen

(liquid) region. The one—dimensional mathematical description of the problem in the frozen

regionis
3’1" 2 3T, 1 3T. 0<r<s(t),t>0 (3.14)
— '0' _ = _
3r2 r37 (1,7
and in the unfrozen region is
3’T. + 2371 g 1 37: s(t)<r<oo.t>0 (3.15)

57 7‘37 as;
where 3(1). the freezing front. is the location of the interface separating the frozen and unfrozen

regions. Thetemperatureisfiniteasr-ioo
7'04) 3 1", r-)°°,t>0 (3.16)

and.forconvenience,thepointheatsinkisassumedtoincreasewiththesquarerootoftime

l9

lim 2 3T, . m r—)0,t>0 (3.17)
r 441" k’??] 2Q(a,t)

The uniform temperature initial condition is

1,03,) 3 1i 0<r<oo,t=0 (3.18)
Diamdnpmsednngeoccmgdunngmefieezingpmcess.mrfacewndifiomueflso
required. A continuity of temperature at the interface requires that

T,(r.t) = T,(r.t) = T. r=s(t). t) O (3.19)

and. from the conservation of energy at the interface

I"??? - 19.335 -.- p L £9 r = 3(1). r> 0 (3.20)
,

To obtain analytical solutions to the describing differential equations. a variable

transformation is introduced. The dimensionless similarity variable 11 is chosen to be

r
= W (3.21)

The location of the freezing front. 3(1). as a dimensionless variable 3.. is assumed to be (Ozisik.

1980)

so)
1. . an
my ( )

Also necessary is the transformation to dimensionless temperature variables as follows:

1,4,.

9 . (3.23)
' n-Tr'
e, = if; (324)

20
The mathematical description of the direct problem. in dimensionless form. is obtained from the

substitution of equations (3.21). (3.23). and (3.24) into equations (3.14) and (3.15). The

dimensionless temperature profile in the frozen region is described by

 

 

:19; +[%+2fl}:%'0 0<r|<2t (3.25)
and in the unfrozen region by
:31 + [7? + mad}? . 0 11>). (3.26)

where a," the dimensionless thermal diffusivity. is defimd to be

cg, . 3.1 (3.27)
at

The dimensionless temperature as n -) co. obtained from the substitution of equation (3.16) into

equation (3.24). leads to the boundary condition
elm—i...) 3 o 11 —i on (3.28)

The second boundary condition is obtained by differentiating equation (3.23) with respect to 11.

and substituting the resulting derivative into equation (3.17). It is expressed as

lim 209’ . . n -) 0 (3.29)
n-i 2““ at] Q
where k... the dimensionless thermal conductivity. is defined as
k“ = 51 (3.30)
*1

and Q‘. the dimensionless heat sink coeflicient. is defined as

21

° = __Q 3.31
Q 2xk,(T.-T..) ( )

The dimensionless interface condition at n = 7t is obtained from equations (3.19). (3.23) and

(3.24). and is given by

e,(n =2.) = em =1) = 1 11 = 1 (3.32)

The dimensionless temperature variables. equations (3.23) and (3.24). are used with equation
(3.20). resulting in the following dimensionless expression for the conservation of energy at the

interface:

09 d9 = it (3.33)
k ,_' - _'] = L a n
[ d" Jud. [d7] not

where L'. the dimensionless latent heat of fusion term. is defined as

u .. .11"— (3.34)
c~(r_-r,.)

The solutions to the final forms of the differential equations describing the temperature

profiles ofthe frozen and unfrozen regions are assumed to be ofthe form (Paterson. 1952)

. 1 «1*- J17 , 0<n<x (3.35)
9,01) A [2'11" Tureen] B
1 «1'3. - (’7 - >7. (3 36
9.01) = C -—-T,,-e —erfc(naa"‘) D '1 ~ )
21101,, 2

The statements for the dimensionless boundary and interface conditions. equations (3.28). (3.29).

(3.32). and (3.33). are used to solve for the four unknown coefficients

22

 

A = _Q . (3.37)
. e 4"" - J; 4 (3.39)

C [23.1132 Tejmfi]

D = 0 (3040)

The final form of the solution for the dimensionless temperature profile for the frozen region is

expressed as

 

 

g - -tl; _ 0 < n < A (3.41)
6,01) 1 Q‘[-ef:‘-1---M ‘2}?- T —(¢IfC(n) arm»)
and for the unfrozen region as
’01.. 41a, '1
6.01) g e a m _ flew. a?) e m - Emmy) n > i (3.42)
2113.1 2 2M“ 2

The derivatives of these dimensionless temperature profile expressions with respect to 11
are used in the dimensionless interface condition. equation (3.33). to obtain the following

transcendental equation for the dimensionless freezing from location:

11 = 1. (3.43)

 

23
resulting in the function M). which equals zero

fl1)'0'k,,Q{%;-]' 11:). (3.44)

 

3.3 Estimation of the Thermal Properties

Estimation of the thermal properties of the tissue to be frozen. using parameter estimation
techniques. requires the calculation of the sensitivity coefficients for each property for use in the
Box-Kanemasu Method. see equations (3.7) and (3.8). A sensitivity coefficient is defined as the
change in a given variable due to a change in a specific parameter, with all other parameters held

constant. Mathematically. it is defined as

X ‘rs- ' l3 [-32-] w (3.45)

where X o”. is the dimensionless sensitivity coefficient. 1" is the process variable (i.e. the
dimensionless temperature 6). B the specified parameter of interest. and E,- ¢ B are all the other
parameters other than B.
3.3.1 Determination of the Sensitivity Coefficients

The thermal properties estimated in this investigation are the dimensionless latent heat of
fusion. L'. the dimensionless thermal conductivity. k... and the dimensionless thermal difiusivity.
01,. The sensitivity coefiicients are obtained by differentiating the dimensionless temperature
profile expressions, equations (3.41) and (3.42). with respect to each parameter (1:. 1:... and (1,).
In the frozen region the sensitivity coefficients are determined fiom the following expressions:

3°: . £90,3—; 0<11<lt (3.46)
av War:
30 39 a), 0 <

. g . n < A (3.47)
it: an:
30 39 a; 0 <

s . a n < A. (3.48)
36: 31'3on

where

39- . '9 ’e 4' (3.49)
3T 21.5

In the unfrozen region the following expressions are obtained:

33; . 39:31; rpx (350)
BL‘ 35531.:
39: g 39:31 11>). (3.51)
37.. 3535

 

1/2

., {[mfl; geese] - [we mean]

 

 

E 41111.3," 231:}? 21101,,
‘C 4'11, C 4%., ax e 4%.: fu- 1 .2 n > x (3.52)
x - - rfCOa ")
Lita? mafia—I; ]i [2241}? TC " i
where
39, .9 ‘1’“: J; m '6 4"" C 4“ __ J; .2 (3.53)
at [ 21px;" + 707.001“: ) ][W m 762000;”)

The partial derivative of l. with respect to each parameter is determined as follows (Scott.
1993): since the tramcendental equation. equation (3.44). is equal to zero. the total derivative of
10.) is also equal to zero

25

3f
. p .1) = 0 =
m g ”a [w] {Mlmw
,2’: [3f] 4‘ , [3f] a (354)
a; a: 3.1m 3x ‘M’m

 

 

Dividing equation (3.54) by dB and solving for deB yields

 

 

 

%__[%J.-.m*§[%]mir
(3%) ...,.........

Neglecting the higher order terms. the partial derivative of A with respect to B can be

approximated from equation (3.55) as

 

ar g _ [3%) rem—- (3.56)
”‘2'“ [3%) ..........

The partial derivative of 10.) with respect to 1. is determined using equation (3.44). and is found

tobe

 

 

2 4s + 4: _ 2 4h, _ 4%. 4%.,
#“*‘Q'[“ )8 . i'H—Menre Hindi” 'éflfmfij

e are, _e an, e are, JI— ‘ '2 (357)
- __ __ __ - _erfc(3.a a) - L '
[Wariizw i} [2103” 2 .. i

26
Differentiating 1'0.) with respect to L', k... and (1,. leads to the following expressions respectively:

 

 

Ll") - -)t (3.58)
31. '
31m , Q 'e ""
'37:": T (359)
- are,
aflx) . - 'flzahne 4%“ " Game 4%., e - fle’f ; Ln)
3an 4120’, 2211:,” T
e are, -e are, e 4%., J; '2 (3.60)
- - —erfc(ka‘”)
[mum] [We 2 "

The sensitivity coefficients required for use in the Box-Kanemasu Interpolation Method for the
estimation of the thermal properties are completely defined by equations (3.46) - (3.53), and (3.57)
- (3.60).

In summary, the partial differential equations governing the temperature profiles of the
frozen and unfiozen regiorn were made dimensionless through the use of a similarity variable
transformation and dimensionless temperature expressions. Analytical solutions were then found
for these dimensionless differential equations. These analytical solutions describe the
dimensionless temperature profiles of the frozen and unfrozen regions of the tissue that is
cryosurgically frozen. From the dimensionless temperature expressions. sensitivity coeflicients
for the thermal properties were obtained for use with the Box-Kanemasu Interpolation Method.

used to estimate these material properties.

27
3.4 The Inverse Problem of Deternrining the Optimal Cryosurgical Treatment Time

The objective in this portion of the investigation was to determine the optimal
cryosurgical treatment time required to achieve a desired minimum temperature at a specified
radius locationofthetumor. Thisdesiredminimum temperaturewilloccuratatimethatis
gremermmmeacmdheaunanfimeduemmedifflnionofanrgymatconfinuesaflerme
treatment has stopped. To account for this continued cooling after the cessation of treatment. the
principle of superposition is used to describe the dimensionless temperature profiles after the

cryosurgical treatment time. 1,. as follows:

em) = em) - om) ‘1 <n' (161)

where 11 is defined in equation (3.21), and

f
. . _ (3.62)
Tl 2(a,(r - :9)“

The Box-Kanemasu Interpolation Method is again utilized for the determination of the
optimal treatment time. Therefore. the analytical solutions to the describing differential equations.
determined in Section 3.2. are once again used. Applying the principle of niperposition. the
dimensionless temperature profile in the fiozen region is expressed as

e ""

a — ' - C 4’ .; fi- —
6,01) 1 Q [—211— ‘27 —2-(erfc(n) erfcml]

.. - - 8",!“ -5 -
1 Q [3?- 7X 7-(erfo(11 ‘) crawl]

0<n<Ln<n° (3.63)

and in the unfrozen region as

28

 

in, 4b, '1
am) - [2“: m - 15.20%] of)] [5251.” - genmfl

as!

_ 94"" _ J; . in in" _ J; IR -1 11>AJI<TI. (3-64)
“2“ m Te'fcm a.) )J [m Tasman )] }

Aswidrdreesfimafionofmemermalpmperties,medeterminanonofmeopfimal
treatment time. t,. also requires the calculation of the sensitivity coefficients for r, from the frozen
and unfrozen region dimensionless temperature profile expressions. The sensitivity coemcient for

the frozen region is obtained by differentiating equation (3.63) with respect to 1,. resulting in the

followingexpression:
39m) 39(n')an' 0< -
s _ . n<Ln<n (3.65)
5!. an' '3?
where
3991') _ -Q'e"‘°a O<n<Ln<n° (3.66)

an ‘ 21)5
Differentiation of the expression for '0‘. equation (3.62). with respect 1. yields

311' r
. (3.67)
7‘: 4af”(r - if”

The sensitivity coefficient for the rmfrozen region is obtained by differentiating equation (3.64)

withrespecttorgitisdeterminedtobe

39m) _ 39.01? an‘ n>Ln<rf (3.68)
T. —8n' "a:

where

39,01.) a e‘“‘ e'm“ _ ft— 1 ‘1 >1, (11° (3.69)
T [W] [—115 Teddmufl) '1 fl

Duemdediifusionofmergymatwndmwsafiertheueannenthassmpped.medenmd
minimumternperatureatthespecifiedlocationwillacmallyoccuratatime.r_,-,.thatisgreater
than the treatment time. r.. Therefore. in the determination of the optimal treatment time it is also
mcessarymcdcmmedeaaualfimemmmespecifiedbcafionmachesmedesimdmmimum
temperature. This is accomplished by minimizing the expressions for the temperature profiles in
the frozen and unfrozen regions and solving for this actual time. t“... i.e.. by setting the partial
derivatives of the temperature expressions with respect to 1 equal to zero. Differentiation of the
dimensionless temperature profile for the fiozen region. equation (3.63). with respect to 1 yields

the following expression:

 

39,00 20:39(11)an 39,'(1‘|')an 0<n<Ln<n° (3.70)
T: m E T an —T
where
39.01) , 9'8 "V 0<n<x (3.71)
317. 2112

and 80,01‘) is given by equation (3.66). Differentiating the expression for ’0, equation (3.21).
an o
with respect to t yields

an . -r
.5}. W (3.72)

while differentiating the expression for n’. equation (3.62). with respect to t results in the

following expression:

3‘" . ___’_'_ (3.73)
7’- 4a}”(r - t)”

For the unfiozen region. the partial derivative of equation (3.64) with respect to 1 yields

39,01) _0 . 30mm + 30(1)) 311° n>Ln<n‘ (3.74)
T: u 5i '3:- atr '3:—
where
aern) _ -e*" a“ _ Jr? in '1 >71 (3.75)
T [WNW T‘”“‘““’] "

and 39,01‘) is given by equation (3.69).
an.

In summary. the Box-Kanemasu Interpolation Method. used in the determination of the
Optimal cryosurgical treatment time. required the determination of the temperature profiles and
sensitivity coefficients for the frozen and unfrozen regions. The dimensionless temperature
expressions used to estimate the thermal properties of the tissue were again used. with the
principle of superposition applied to describe the dimensionless temperature profiles of the tissue
after the cryosurgical treatment time. The required sensitivity coefiicients for t, for both the
frozen and unfrozen regions were calculated Also necessary was the calculation ofthe actual
time. r...mwhichdredesimdmmimumtanperammisachievedatmespecifiedlocafionduem

mediffusionofenergydratcontimresafiermecessafionoftheueaunem.

CHAPTER 4

ANALYTICAL PROCEDURE

Inthischaptertheproceduresusedtoestimatethethermalpropertiesofthetissuetobe
cryosurgically frozen and to determine the optimal treatment time required to achieve a desired
mhimum¢mperatumflaspedfiedbcafionuedescfibed1hresdtsobminedfiommis

investigation are presented in Chapter 5.

4.1 The Sensitivity Coefficient Analysis for the Thermal Properties

Prior to the use of the Box-Kanemasu Interpolation Method to estimate the thermal
propertiesofthetissuetobedestroyedbycryosurgicallyfreezingitwasimportanttoexaminethe
sensitivity coefficients for magnitude and linear dependence. A sensitivity coefficient with a
-all magnitude (<10’) indicates that the dimensionless temperature profile is relatively
insensitive to changes in a given parameter. while a large magnitude (21) indicates extreme
sensitivity to a change in a specified parameter. A sensitivity coefiicient with a small magm'tude
indicates that there is very little information about the value of the parameter available hour the
temperature measurement data. making estimation of that parameter dificult or impossible.

Another important consideration when using the Box-Kanemasu Interpolation Method is
the possibility of correlation existing between the parameters to be estimated. To simultaneously

estimate two or more parameters. it is necessary that their respective sensitivity coefficients not

31

32
be linearly dependent. If the sensitivity coefficients are found to be linearly dependent. the
parameters are correlated and cannot be estimated simultaneously.

To conduct the analysis of the sensitivity coeflicients for the thermal properties to be
estimated. a Fortran program. SENSEFOR. was written. The sensitivity coefficients were
determined using equations (3.46) - (3.53) and (3.57) - (3.60) for each property under
consideration. To calculate the sensitivity coefficients. it was necessary to assign values to the
dimensionless heat sink coefficient and thermal properties. The dimensionless heat sink
coefficient was kept constant at a value of

Q‘ = -1.0
The thermal properties to be estimated. the dimensionless latent heat of fusion. L'. the
dimensionless thermal conductivity. k”, and the dimensionless thermal diffusivity. or,,, were
assigned values of:

L' -100.0

k“ = 1.0

(1,, = 1.0
which remained unchanged throughout the entire investigation.

The sensitivity coefficients for the thermal properties were calculated as functiorn of the
independent variable 11. The independent variable was varied over a range of 0.01 to 2.0 in steps
of 0.01. For the above values of the dimensionless heat sink coefficient and thermal properties.
the location of the dimensionless freezing front 3. was determined by solving the transcendental
equation. equation (3.44). and found to be 0.17302 (using the root-finding subroutine ZBRENT‘.
Press et al.. 1986). Therefore. the range of 1| sufficiently covered both the frozen and unfrozen
regions. The program provided an output of n and the sensitivity coefficients for the thermal

properties in a nondimensional form. expressed as

33

Xu . L ._ (4.1)
3L '
w J
c ‘ (4.2)
x,‘ In}:
Xa‘ s (1% (4.3)
I

These dimemionless forms of the sensitivity coefficients are plotted versus 11 in Figure 4.1. A

copyofthisprogramandasampleoutputfilemaybefoundinAppendixA.

.0
or

 

4—..‘

.0
U
9;

.0
.a
J
I
I

 

 

 

 

 

-o.3- i
: — x
l a.
.mi " x»
--- It.
-Oo5 I I 1 T r r 1 I 1 r r 1 r v fir Iij I r I r
0.0 0.4 0.8 1.2 1 .6 2.0

Dimensionless Sensitivity Coefficients

77

Figure 4.1. Dimensionless Sensitivity Coemcients versus 11

4.1.1 Magnitudes of the Sendtivity Coefficients
Figure 4.1 demonstrates the magnitudes of the sensitivity coefficients in the frozen and

unfrozenregionsandhowtheychangeattheinterface. Themagninldesofthedimensionless

34

sensitivity coefficients XL, and x“ are the largest within the fiozen region. They decrease at
theinterfacebenveendrefrozenandtmfiozenregionsandapproachzeroasnincreaseswithin
the unfrozen region. This indicates that the most available information about the values of L' and
kdiscontainedinthesimulatedtemperanrre measurementdataobtainedfromthefrozenregion.
with very little information available from the unfiozen legion data

Incontrast.themagnitudeofthedimensionless sensitivitycoefficient xa‘ isquitesrnall
in the frozen region. It changes sign at the interface and increases slightly in the unfrozen region
before approaching zero as 11 increases. Therefore. the most available information about the value
adiscontainedinthesimulatedtemperamremeasurememdataobtainedfiommeporfionofme
unfrozen region adjacent to the interface.

0f the three dimensionless sensitivity coefficients plotted in Figure 4.1. XL. has the
largest magnitude. followed by KL and x“; This indicates that the temperature measurement
data provides more information about the value of L' than it does of k,‘ and 0... Therefore. of the
three thermal properties being estimated. the estimate of L‘ should be the most accurate. followed
by k, and a“.
4.1.2 Linear Dependence Between Sensitivity Coefficients

The issue of correlation existing between parameters was addressed by plotting one
dimensionless sensitivity coefficient versus another to observe any linearity between them. The
following graphs were generated: X,‘ versus XL” x9. versus XL” and xc‘ versus ix“,
These graphs are presented in Figures 4.2. 4.3. and 4.4 respectively.

The thermal properties L‘ and I:.' are correlated throughout the frozen and unfrozen regions
of the tissue. as demonstrated by the linear relationship between their sensitivity coefficients in
Figure 4.2. Thedashedlineinthisfigurerepresents adiscontinuityattheinterface. Therefore.

these parameters could not be estimated simultaneously using the Box-Kanemasu Method.

35

 

 

 

 

 

 

 

 

‘5
.2 0.5
O
afi ‘
“a
o 0.44 ‘~.
0 ~.\
A ‘ ‘s‘
,3 ~.
as 0.3“ ‘\‘
f: n.
m ‘s
.e' .
0 s
m 0.2- “\
m
In .
o
"a
2 0.1-
m
a ~1
a)
,§ 0.0 . . . .
Q -O.5 -O.4 -0.3 -O.2 -0.1 0.0
Dimensionless Sensitivity Coefficient
Figure 4.2. Dimensionless Sensitivity Coefficiernts xi versus XL‘
‘5
.2 0.04
‘8 0.02.l ‘s‘
c . ‘~.
0 x‘
g‘ o.oo- x‘
:*-’ . ‘~.
.E
3:: -0.02- ‘ ‘4
m 8
8x3 .
m -0004-
g d
a) -0.06-
III.
a «l
c
'3 -O.08d
‘3 .
a)
.g -o.lo i I 1 I T r I r
g -o.5 -o.4 -o.3 -o.2 —o.1 0.0

Dimensionless Sensitivity Coefficient

0

Figure 4.3. Dimensionless Semitivity Coefficients X“. versus xv

36

The thermal propertiesL'andadare correlated withintlnefrozen region ofthe tissue. as

shownbythesinglepointinFigure4.3. Againtlnereisadiscontinuityattheinterface.

represerrtedbythedaslnedline. Intlreunfrozenregionthethermalpropertiesarenotcorrelated.

asindicatedbythenlrvedportionofthisfigure. However.asnincreasesandthesemitivity

coefficieras approach zero. the properties becomes correlated again.

y Coefficient

X

Dimensionless Sensitivit

Figure 4.4. Dirnernsionless Sensitivity Coeficieras 1‘. versus x1Ir

0.04

 

0.02-
0.00-

-0.02~

‘d

-0.04-
-0.06-r

-0.08-

 

-0.10

 

 

0.0

r
0.1

Dimensionless Sensitivity Coefficient

r
0.2

I
0.3

Ti

1
0.4

0.5

Figure4.4denonstratesthattheparametersk,,andn,,exhibitthesamebehaviorasL'and

a... therefore. it was determined tlnat L‘ and a. could be estimated simultaneously. as could It.

myprovidedmemfiiciemtenpeammmeasuremendnafiommemfiozmmgimwasmed.

37
4.2 The Sum of Squares Function for the Thermal Properties

Prior to the estimation of the thermal properties using the Box-Kanemasu Interpolation
Method.itwasofinteresttoexaminethesmnofsquaresftmctionforeachofthesethermal
properties. ShrcemesumofsquaresfuncfionismbemhfimizedinthismeMaflatminimum
would indicatethattheparticularparameterwouldbemore difficulttoestimatethanasteep
minimum, requiring more iterations for convergence. A Fortran program, SQUARESPOR, was
written using the dimensionless temperature profile expressions for the frozen and lmfrozen
regions, equations (3.41) and (3.42) respectively. A copy of this program is provided in Appendix
B. Using the assigned values for the thermal properties this program calculated the dimensionless
temperature values as r] was varied from 0.01 to 1.5 in steps of 0.01. These values represent the
temperature measurement data vector Y in the sum of squares function, equation (3.1). The value
of one thermal property was then varied in -all increments, with the dimensionless temperature
values calculated as n was again varied from 0.01 to 1.5. These values represent the calculated
temperamredatavectorTinthesumofsquares function. 'Ihesumofsquaresfunctionwasthen
calculated for each increment of the thermal property. and is plotted versus the varying thermal
property for L' ill Figure 4.5, and for ltd and adin Figure 4.6.

Asdemonstratedbythesefigures,thesumofsquales functionforadhasthesteepest
minimum. This indicates that the estimation of a... should require the fewest number of iterations
for convergence in the Box-Kanemasu Interpolation Method. The minima of the sum of squares
functions for L‘ and ks are considerably less steep, indicating that the estimation of these thermal

properties would require more iterations for convergence.

38

 

5.0E-07

q

4.0E-07-

q

JOE-0'7-i

2.0E-07q

1.0E‘07“ ‘\ I

'4 ‘ I

Sum of Squares

 

0.0E+00

1

 

 

-1.0E-07 .
-160

 

I ' I I 1
-140 -120 -100 -80

Dimensionless Latent Heat of Fusion L’

Figure 4.5. The Sum of Squares Flmction versus L‘

I
-60

-4o

 

5.0E-07

q

4.0E-07d ‘.

"‘ s

3.05.07“ “s.

2.05-07-
4 \‘
1 .05-07‘

J
\\. \ I ,4"

Sum of Squares

 

0.0E+00

q

 

-1 .OE-07 1

0.4

I
0.6

I
0.8

I
1.0

I
1.2

Pigure4.6. TheSllrnoquuar'erll'trrlctionversusk,,andtxu

31"

1.4

 

1.6

Dimensionless Thermal Conductivity k.
Dimensionless Thermal Diffusrvlty a.I

39
4.3 Estimation of the Thermal Properties

A Fortran program of the Box-Kanemasu Interpolation Method was used for the
estimation of the thermal properties. This program. NLINAFOR (Beck. 1991), required
modification prior to use. The first modification involved the use of the dimensionless
temperature profiles for the frozen and unfrozen regions, equations (3.41) and (3.42) respectively.
in the subroutine MODEL. These expessions provide the calculated temperature data necessary
for use in the Box-Kanemasu Interpolation Method. as required in equation (3.8).

The second modification involved the placement of the expressions for the sensitivity
coefficients for the thermal properties. equations (3.46) - (3.53) and (3.57) - (3.60). into the
subroutine SENSE. This subroutine calculates the sensitivity coefficients for each parameter under
consideration to be used in the Box-Kanemasu method for the estimation of the thermal properties
as required in equations (3.7). (3.8), and (3.11). A capy of this program can be found in
Appendix C.

4.3.1 Input for the Box-Kanemasu Interpolation Method

The use of the Box-Kanemasu Interpolation Method for the estimation of the thermal
properties required an input of internal temperature measurement data. To simulate this data, a
Fortran program was written. This program, MODFOR. uses the dimensionless temperature
profile expressions. equations (3.41) and (3.42), and the assigned values of the dimensionless heat
sink coefficient and thermal properties to calculate the values of temperature as a function of the
independent variable 1]. with n varying from 0.01 to 1.5 in steps of 0.01. These calculated
temperature values were used to simulate the internal temperature measurement data from the
frozen and unfrozen regions as required in equation (3.8). Also included in this program was the
subroutine RANDOM (Press et al., 1986). a random number generator used to simulate random

measurementenors. Auserinputofthestandarddeviationoftherandomnumbersisrequired.

40
thereby allowing the standard deviation of the measurement errors to be varied. Another required

user input is the seed or initialization number, each different negative number produces a different
setofrandomnumbers. Thesesimulatedmeasurementerlorscouldbeaddedtothesimulated
measurementdata. AcopyofthisprograrnandasampleoutputfilearefoundinAppendix D.
4.3.2 Individual Estimation of the Thermal Properties

Initially. the estimation of the thermal properties was conducted on an individual basis.
Using the simulated temperature measurement data provided by MODFOR as input for the
modified version of NLINA.FOR, the first property to be estimated was the dimensionless latent
heat of fusion. L’. This property was estimated using exact temperature measurement data. i.e..
without measurement errors. To estimate L‘ with measurement data that contained errors three
measurement error standard deviations were used: 0.1. 1.0 and 10.0. For each standard deviation,
twelve different sets of random measurement errors were added to the simulated temperature
measurement data and twelve estimations of L' were conducted. Since the actual value of L‘ was
-100.0. an initial estimate of -50.0 was used for the first six estimations; and an initial estimate
of -150.0 was used for the remaining six estimations.

To include prior information of the value of L’. the sum of squares frmction was modified
as in equation (3.4). This required slight revision of the subroutines MODEL and SENSE in the
NLINAFOR program. Copies of these subroutines may be found in Appendix C. The input file
of simulated internal temperature measurement data was also modified slightly, a c0py of this file
is located ill Appendix D. Using a standard deviation of 0.1 forthe prior information, L‘ was
estimated using exact temperature measurement data. Twelve estimations were performed at each
of the three measurement error standard deviatiom. Initial estimates of -50.0 and 450.0 were
again used. This approach was repeamd using prior information standard deviations of 1.0 and
10.0.

Themethodusedfortheindividual estimationofthedimensionlessthermal conductivity.

41
lg, and the dimensionless thermal diffusivity. or“. was the same as the one described for L’ with

the following exceptions: I) the standard deviations used for the measurement errors were 0.001.
0.01. and 0.1: 2) the prior information standard deviations were also 0.001, 0.01, and 0.1; 3)
because the actual values ofbothlt,I and or" were 1.0. initial estimates of 0.5 and 1.5 were used.
4.3.3 Simultaneous Estimation of the Thermal Properties

As demonstrated in Figure 4.2. the thermal properties L‘ and It, are conelated. thereby
eliminating the possibility of simultaneous estimation of these two parameters. However. as
shown in Figure 4.3. the parameters L' and agate uncorrelated in the unfrozen region and could
be estimated simultaneously. From Figure 4.4, it was determined that lg, and cg, are not correlated
in the unfrozen region and could also be simultaneously estimated.

The simultaneous estimation of L‘ and a“ began without the use of prior information. The
properties were estimated using exact temperature measurement data. To estimate the thermal
properties using measurement data that contained errors, three standard deviations of measurement
errors were used: 0.0)]. 0.01 and 0.1. For each standard deviation. twelve different sets of
random measurement errors were added to the simulated temperature measurement data and twelve
simultaneous estimations of L‘ and a... were conducted. Initial estimates for L' and (1,, of -50.0
and 0.5 respectively were used for the first six estimations, and -150.0 and 1.5 respectively were
used for the remaining six estimations.

To include prior information of the values of this pair of thermal properties, the modified
sum of squares function. equation (3.4), was again used. Using prior information standard
deviations of 0.1 for L‘ and 0.001 for a... the thermal properties were estimated using exact
temperature measurement data. Twelve estimations were then performed at each measurement
error standard deviation: 0.001, 0.01, and 0.1. Again. initial estimates of -50.0 and 0.5 were used
for the first six estimations. and -150.0 and 1.5 were used for the remaining six estimations. This

approach was repeated using prior information standard deviations of 1.0 and 10.0 for L‘. 0.01 and

42

0.1 for a“.

The simultaneous estimation of ha and (1,, followed the same procedure as L‘ and or...
with the following exceptions: 1) prior information standard deviations were 0.001, 0.01. and 0.1:
2) because the actual values ofk,, and a, were 1.0, initial estimates of 0.5 and 1.5 were used.

For each set of twelve estimations, the estimated thermal property values were averaged.
with the standard deviation and a 95% confidence interval calculated. A comparison could then
be made between the estimated thermal property values provided by NLINAJ-‘OR and the actual
property values used to generate the simulated temperature measurement data. The number of
iterations required for convergence was also averaged for each set of twelve estimations. with the

standard deviation calculated.

4.4 The Determination Procedure for the Optimal Cryosurgical Treatment Time

To destroy undesirable biological tissue by cryosurgically freezing it is of extreme
importance to be able to accurately determine the optimal cryosurgical treatment time required to
provide a desired minimum temperature at a specified location. not only to ensure destruction of
thediseased tissue.butalsotominimizelossoffllesurroundinghea1thytissue. lnthissection.
theproceduleusedtodetermine thisoptimaltreatrnenttimeispresented.

4.4.1 Sensitivity Coefficient Analysis for the Optimal Cryosurgical Treatment Time

Because the optimal treatment time was to be determined individually, linear dependence
between sensitivity coefficients was not a concern in this portion of the investigation. However.
it was of interest to examine the magnitude of the treaunent time sensitivity coeflicients
determined from the frozen and unfrozen regions. To accomplish this, the program SENSEFOR
was modified by replacing the thermal property sensitivity coefficients with the expressions

formulated in Chapter 3 for the treatment time sensitivity coefi'rcients. equations (3.65) - (3.69).

43

Thetreatmenttimewasassignedavalueof

rc = 0.1850 seconds
Although the value of the treatment time may not be a reasonable duration of treaunent. the
objective of this investigation was to assess the minimization procedure rather than to simulate
an actual cryosurgical procedure.

Theacmalfimeatwhichdledesimdmmumtanpemmreisachievedataspecified
location due to the diffusion of energy after the cessation of treatment (determined from the
program MODC.FOR. discussed below) was

r... = 0.1865 seconds
Thiswndnuedwofingofflnfismafierdnardofflwayomrgicflfieammtisdemonmw
in Figure 4.7. A dimensionless temperature greater than 1 indicates that the tissue is frozen. while
aternperanlrelessthanlmeansthatthetissueisunfrozen.

To calculate the sensitivity coefficients for re. the independent variables 11 and n‘ were
varied by varying the value of the radius from 0.01 to 0.5 in steps of 0.01. The program provided

an output of n and the dimensionless form of the sensitivity coefficients for re. expreswd as
x .. t 3°" (4.4)

To observe the magnitude, the dimensionless sensitivity coefficient values for t. were
plottedversusn:thisgrafllispresentedinFigure4.8. Asdemonstratedbythisfigure,the
magnitude of the dimensionless sernitivity coefficient for the optimal treatment
time. X”, becomes very large as it approaches zero. Therefore. the most available information
ofthe value of r. is contained in the dimensionless temperature measurement data obtained from

the frozen region.

 

2.5

‘ Minimum Temierature Achieved at t...

2.41

1

2.3-1

   

“I... = 0.1865

2.2e

Dimensionless Temperature

‘ t. - 0.105»

 

 

I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
1

 

 

2.1.,.,..j...,.,.T.l
0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25

Time (seconds)

Figure 4.7. Dimensionless Temperature of the G'yosurgically Frozen Tissue versus Time
at a Radius of0.l Meters. Given aTreatment Time 2, = 0.185 Seconds

 

 

 

 

 

 

.9
t:
.2 250
O
2:: « 1
"o‘
o 200‘
U 1
3 1 50
IS
:0: 1
§>g 100‘
U) ..
3 50~
0
F!
g d
ea 0.:
a d
0
.g ‘50 ' I ' I ' I ' I ' I '
C: 0.0 0.1 0.2 0.3 0.4 0.5 0.6

Figure 4.8. Dimensionless Sensitivity Coeficierl 1:. versus 11

45
4.4.2 The Sum of Squares Function for the Optimal Cryosurgical Treatment Time

Pliortothedeterminationoftheoptimaltreatmenttime,itwasagainofinterestto
examhndnmofsquaresfiulcfimfortheneamemfime.r,mdaemineifdnminhnumof
thisfunction was fiatorsteep. TheprogramSQUARESFORwasagainuaedafterslight
modificafimbyusingmeexpresdomfmmedimmsioMesswmpaaumpromesfordnfiown
and unfiozen regions after the treatment time, given by equations (3.63) and (3.64) respectively.
Themmofsquamsfiuwdonwascakuhtedfolbwmgmenmeprocedundesuibedeecfim
4.2withrcvariedinsmallincrements. Thesumofsquaresfunctionisplottedversustfinfigule
4.9. Thisfigure demomtratesthatthesumofsquares flmctionhasaverysteepminimum.
especifllywhenmmpmedmdumofsquuesfimfiomforflndmmflpropemes,mdicafing
that convergence in the Box-Kanemasu Interpolation Method should require a mall mlmber of

iterations.

5.00E-03

 

q

4.00E-03-e

I
I
I
I
I
I
I
I
‘.
3.005-03~ ',
‘.
‘I
2.005-03- 1‘
I
I

1.005-03- \ 'l

‘ I

0.00E+00 3';

Sum of Squares

 

 

 

0.12 0.14

 

r I m m ' I ' I fl I '
0.16 0.18 0.20 0.22 0.24 0.26
Treatment Time t.

Figure4.9. ‘I‘heSmnoquuaresFuncticnversust,

46

4.4.3 Modification of the Box-Kanernasu Interpolation Method Program

TheFomanpmgmmNLmAFORwasalsousedmdeterminemeopdmalueannaufime
with the following modifications: the dimensionless temperature profiles for the frozen and
unfrozen regions after the treatment time, equations (3.63) and (3.64) respectively, were used in
the subroutine MODEL. In the SENSE subroutine, the sensitivity coefficient expressions for r,,
equations (3.65) - (3.69). were used. Samples of these subroutines are located in Appendix E.
4.4.4 Input for the Box-Kanemasu Interpolation Method

Using a cryosurgical treatment time of r‘ = 0.185 seconds. it was necessary to determine
the radius locations at which the desired minimum temperatures were achieved to be used as input
for NLINA.FOR. At the boundary of the tumor the desired minimum temperature was chosen to
be 2.39414, this corresponds to a treatment temperature of approximately -50°C and an initial
temperature of 37°C. Since it is not desired to freeze the sunounding healthy tissue. a second
desired minimum temperature of 0.99467 (approximately 0.15°C) was chosen to be achieved at
amalldistancebeyondthe tumorboundary. Todeterminetheladius locationsatwhichthese
desired minimum temperature were achieved the program MODCFOR was written. In this
programtimewasvariedinsrnall increments. Whentimewaslessthanthetreatmenttime
equations (3.41) and (3.42) were used to describe the temperature profiles in the frozen and
tmfrozen regions respectively. When time was greater than the treatment time equations (3.63)
and(3.64)wereusedtodescribethetemperature profiles. Alsovariedinthisprograrn.insmall
increments. was the radius. The output of this program was dimensionless temperature values at
corresponding radii for given times. This allowed for the determination of the actual time the
minimum temperatureoccursatagivenradiustobet...=0.l865secondsforatreatmenttime
of 0.1850 seconds. From this output, the desired minimum dimensionless temperature of 2.39414

wasachievedataradiusofOJOOmeters. 'I‘hisrepresentstheboundaryofthenlmor. Thesecond

47
desired minimum temperature of 0.99467 was determined to be achieved at a radius location of

0.1500 meters. Although the size ofthe radii are large, due to the chosen values ofthe thermal
properties, it was not of importance since the objective of this investigation was to test the
estimationprocedure forreIianceand accuracyratherthantopredictthetreatmenttimeforan
acmalsituation. AcopyofthisprogramandasampleoutputfilearelocatedinAppendixF.

lnthedeterminationoftheoptimalneatmenttimetheYvectorinthesumofsquares
flmction. equation (3.1), contains the desired minimum temperatures to be achieved at specified
radius locations rather than measured temperature values. The measurement data is now
consideredtobethelocationorradiusofthetissueatwhichthedesiredminimumtemperature
is to be achieved. In the medical setting this information is most commonly obtained by
ultrasonic measurement. and therefore contains measurement errors. Using the desired minimum
dimensionless temperatures and corresponding radius values determined above. an input file of
exact measurement data was generated for use in the NLINAFOR program. A sample of this
input file, TEMP.DAT. is located in Appendix G.

To add random errors to the simulated measurement data. i.e.. the radii. the Fortran
program RADFORwas written. Thisprogram wasdesignedtoreadtheinputfileofexactdata
and, again using the random number generator RANDOM, produce an output file with -ulated
random measurementerrors addedtotheradius values. Thisfile was alsousedasinput forthe
NLlNAFORprogram. CopiesofthisprograrnandouqartfilearealsolminAppendixG.
4.4.5 Determination of the Optimal Cryosurgical Treatment Time

The objective of this portion of the investigation was to determine the optimal treatment
timemmndwachieveadesimdminimumtempemmmatagivmmdiusofdiseasedfismem
bedestroyedbyfreezing. whileresultingintheleastpossibleamountofdamagetothe
sunoundinghealthytissue. Usingdataobtainedflomboththefrozenandunfrozenregionsas

irlputtlleprogramNLmAmeasusedtodetermirletheoptimal treatmenttimenc. 'Ihedesired

43
minimum dimensionless temperature for the frozen region was 2.39414 at a specified radius

location of 0.100 meters. The desired unfrozen region temperature was 0.99467 at a specified
radius location 0.1500 meters. Without the use of prior information 1. was first determined using
exact data. i.e.. without measurement errors added to the radius locations. To estimate rc using
radius measurement data that contained errors, three measurement error standard deviations were
used: 0.0001. 0.001. and 0.01. For each standard deviation. twelve different sets of random
measurement errors were added to the radius measurement data and twelve estimations of r, were
conducted. Since the actual value of rewas 0.185 seconds. an initial estimate of 0.125 was used
forthefirstsixestimations:aninitialestimateof0.2.45 wasusedfortheremainingsix
estimations.

To include the use of prior information of the treatment time. obtained from a previously
performed procedure with the same tumor radius. the sum of squares flmctiorl was modified as
in equation (3.4). Using prior information with a standard deviation of 0.0001, the determination
of r, was performed using exact data Twelve estimations were conducted at each radius
measurement error standard deviation of 0.0001. 0.001. and 0.01, with initial estimates of 0.125
and 0.245 seconds again used. This approach was repeated using prior information standard
deviations of 0.001 and 0.01.

From the sensitivity coefiicient analysis. it was determined that data obtained from the
frozen regioncontainedthemostinforrnationabouttheactualvalue oft... while dataobtainedfrom
the unfrozen region provided less. Therefore. determination of t, was performed using input data
obtainedentirelyfromthefrozenregiontoseeofthe accuracyofthe estimates couldbeimproved.
The program MODC.FOR was again used to determine corresponding desired minimum
temperatures and radius locations within the frozen region. The desired minimum dimensionless
temperatures were chosen to be 2.00999 and 1.41992. corresponding to radius locations of 0.1100

and 0.1300 meters respectively. The estimation procedure. both with and without prior

49
information. was repeated. To determine how accurately rc could be estimated from lmfrozen

region data. the estimation procedure was repeated using input data obtaimd entirely within the
unfrozen region. The desired minimum dimensionless temperatures were chosen to be 0.99467
and 0.83825. corresponding to radius locations of 0.1500 and 0.1700 meters respectively.

4.4.6 Use of Prior Information Obtained from a Different Radius

In the determination of the optimal treatment time. a more feasible source of prior
information is obtained from a previously performed prowdure with a different tumor size. In
the previous procedure. the desired minimum temperature is considered to be unchanged. but the
radius location and the treatment time necessary to achieve that temperature are different. Using
a treatment time of 0.165 seconds the desired minimum dimensionless temperature of 2.39414 was
determined to be achieved at a radius location of 0.0944 meters. To incorporate this prior
information into the procedure used to determine the optimal treatment time, several modifications
of the subroutines MODEL and SENSE of the NLINAFOR program and input file were required.
AcOpy ofthese modified subroutinesandintxrtfilemaybefoundinAppendixH.

Using the prior information obtained from a different radius the treatment time required
to provide a minimum dimensionless temperature of 2.39414 at 0.100 meters and 0.99467 at
0.1500 meters was determined. Using a prior information standard deviation of 0.0001. rc was
firstdeterminedusingexactradiusmeasurementdata. Radiusmeasurementerrorswithstandald
deviations of 0.0001 . 0.001. and 0.01 were used. For each measurement error standard deviation.
twelve different sets of random measurement enors were added to the radius values. with twelve
estimations of I. conducted. Initial estimates of 0.125 and 0.245 were again used for t, This
approach was repeated using prior information standard deviations of 0.001 and 0.01.

4.4.7 Use of Prior Information Obtained from Two Different Radii
Information obtaimd from two previously performed procedures was used to determine

iftheaccuracyofthetreatmenttimeestimatescouldbeimproved. Thedesiredminimum

50
temperature remained unchanged at 2.39414. The first prior treatment time was chosen to be

0.165 seconds. The minimum temperature was achieved at a radius of 0.0944 meters. The second
prior treatment time was chosen to be 0.225 seconds with the minimum temperature achieved at
a radius of 0.1103 meters.

Usingmepnormformafionobtainedfiomthetwodifl’erunraditmeneaunemfime
required to provide a minimum dimensionless temperature of 2. 39414 at 0.100 meters and 0.99467
at 0.1500 meters was again determined. Using a prior information standard deviation of 00101
t, was first determined using exact radius measurement data. Random measluement errors were
then added to the radius values using standard deviations of 0.0001, 0.001, and 0.01, with twelve
estimations of r. conducted at each measurement error standard deviation. Initial estimates of
0.125 and 0.245 were again used for t, This approach was repeated using prior information
standard deviations of 0.001 and 0.01.

For each set of twelve runs. the determined values of the optimal cryosurgical treatment
time. re. were averaged, with the standard deviation and a 95 % confidence interval calculated. The
number of iterations required for convergence was also averaged with the standard deviation
calculated. Comparison could then be made between the optimal treatment time determined using
NLINAFOR and the actual value of the treatment time used to generate the simulated

measurement data.

CHAPTER 5

RESULTS AND DISCUSSION

In this chapter the results obtained for the estimated thermal properties and optimal

treatment time are presented and discussed. The conclusions drawn from these results are

presented in Chapter 6.

5.1 Estimation of the Thermal Properties

In the first portion of this investigation. the dimensionless latent heat of fusion, L‘. the
dimensionless thermal conductivity, k“, and the dimensionless thermal diffusivity. or“. were
estimated both individually and simultaneously.

5.1.1 Individual Estimation of the Thermal Properties

Without the use of prior information. the thermal properties L', k“, and (1,, were estimated
using exact temperature measurement data obtained from both the frozen and unfrozen regions and
with measurement data containing random errors. Prior information of the actual values of these
parameters was then included in the estimation procedure. and the properties were again estimated
using exact data and data containing measurement errors. Three different standard deviations for
both the measurement errors and the prior information were used. The results are presented in

Tables 5.1, 5.2. and 5.3 for L'. k,” and (1,, respectively.

51

52

Table 5.1. Estimation of the Dimensionless Latent Heat of Fusion. L'

Stmdard Deviation ofMeasllemeut Errors

 

0.1

L. I 89.9”
t 0.014
‘buror . 0.010

1.0

L. I JULW
:1: 0.184
fierror - 0.024

10.0

L. 8 ~1m.812
:1: 1.075
%errcr . 0.812

 

L' 8 -99.993

1 0.012
‘baror a 0.1117

L. 8 400.001
5: 0.013
‘Ierror 8 01”]

L. 8 400.001
:1: 0N1
‘berror I: 0.001

 

L' a- -99.990
:1: 0.014
‘baror a 0.010

U s 400.019
:1: 0.160
‘lzerror - 0.019

U - -100.058
:1: 0.078
%error a 0.058

 

 

L’ a -99.990
:1: 0.014
%error a 0.010

 

L. I 400.023
3 0.183
iterror a 0.023

 

L‘ . -100.713
1 0.948
“terror = 0.713

 

Table 5.2. Estimation of the Dimensionless Thermal Conductivity. It,

 

kd 3 LIMI
:1: 0.(X)126
%aror = 0.021

r, .. 0.99782
s 0.01541
m = 0.218

 

k‘ ' LMI
1 0.00008
$610! 8 0.11)]

k, :- 0.99999
1 0.001111
%error = 0.001

 

k, - 1.00017
:1: 0.111108
‘berror :- 0.017

k‘ 3 0.99983
:1: 0.00100
‘baror a 0.017

 

 

 

k. :- 1.00022
1 0.00124
‘laror = 0.022

 

k. 8 0.99799
:1: 0.01345
m = 0le

 

 

 

53

Table 5.3. Estimation of the Dimensionless Thermal Diffusivity. a,

Standard Deviation of Measurement Errors

 

0.(X)I

a, I MIXES
t 0.00115
m = 0.035

0.01

a, - 0.99920
t 0.01120
%error a 0.080

0.1

a, = 1.06361
1 0.05186
%errcr = 6.361

 

ads LW

1 0.00027
%error = 0.1118

(1,‘ I 1.00000
:1: 0.001113
%error = 0.010

a, - 1.00010
1: 0.001110
%error a 0.(X)0

 

a, =- LW
:1: 0.00112
%error a 0.034

a, = 0.99976
:1: 0.00267
%error = 0.024

:1: 0.00014
%aror = 0.017

 

(1,, 2 1.00035
:1: 0.00116
%error = 0.035

(1,, 2 0.99920
1 0.01085
%error a: 0.080

a, 8 1.01299
:1: 0.01085
%errcr a 1.299

 

 

 

 

 

As shown inTable 5.1, estimationofL' using exact measuremerltdatawithoutthe use of
prior information resulted in a highly accurate estimate containing only 0.001% enor, as expected.
With the addition of random measurement errors with standard deviations of 0.1. 1.0. and 10.0.
therewasanoverall decreaseintheaccumcyoftheesfimates.withanassociatedhueaseindle
corresponding 95% confidence intervals. The maximum amount of enor was contained in the
estimate obtained using measurement errors with a standard deviation. 0. of 10.0 and was
determined to be 0.812%. With the use of prior information with a standard deviation of 0.1. the
inaccuracy ofthe estimates decreased, as did the corresponding 95% confidence intervals. The
estimate obtained using measurement errors with a standard deviation of 10.0 contained only
0.001% enor. As the standard deviation of the prior information was increased to 1.0 and 10.0
therewasanoveralldecreaseintheaccuracyoftheestimatedvaluesofL'. Withastandard
deviation of 10.0 for both the prior information and measurement errors. the estimate contained

0.713% error, slightly less than the value obtained without the use of prior information.

54

Without prior information. the estimation of k“ with exact data provided a very accurate
estimate. containing only 0.006% error. as shown in Table 5.2. With the addition of random
measurement enors with standard deviations of 0.001. 0.01, and 0.1 there was an overall increase
inboththeinaccuracyoftheestimatesandthe95% confidenceintervals. Theestimateobtaincd
using measurement enors with a standard deviation of 0.1 provided an estimate of k“ containing
0.228% enor. The use of prior information with a standard deviation of 0.001 resulted in an
hwreasemduaccumcyofmeesfimateswimadecnasemflnassociated95%wnfideme
intervals. As the standard deviation of the prior information was increased to 0.01 and 0.1. there
was an overall decrease in the accuracy of the estimates. When measurement errors and prior
information with standard deviations of 0.1 were used, the estimate was slightly better than the
one obtained without prior information. containing 0.201% error.

As demonstrated in Table 5.3, the estimation of 0,, followed the same trend as the
previous parameters. The estimated value obtained with exact data and without prior information
was highly accurate. containing only 0.002% enor. as expected. The inclusion of random
measurement enors with standard deviations of 0.001, 0.01. and 0.1 in the estimation procedure
resultedinanoveralldecreaseintheaccuracyoftheestimates. Thevalue obtainedusing
measurement enors with a standard deviation of 0.1 provided a very inaccurate estimate containing
6.361% error. The use of prior information with a standard deviation of 0.001 provided a decrease
in the amount of error present in the estimated values of or,” especially when large measurement
enors (o = 0.1) were present. As the standard deviation of the prior information increased. the
accuracyoftheestimatesdecreased. Theuseofmeasurementerrorsandpriorinformationwith
standard deviations of 0.1 provided an estimate of or. comaining only 1.299% enor, considerably
less than the value obtained without the use of prior information.

As expected from the sensitivity coefficient analysis, the estimate of L‘ was the most

accuratewhenexactdatawasusedwithoutpriorinformation. Inthiscasetheestimateofog,

55
contained the next level of accuracy. followed by k, When large measurement errors were

present. the estimate of I, was the most accurate. followed by L‘. The estimate obtained for a.
when large measurement errors (a = 0.1) were present was highly inaccurate. The inclusion of
prior information in the estimation procedure provided an overall improvement in the accuracy of
the measurements. When the standard deviations of the prior information and measurement enors
were large (a = 10.0 for L‘. 0.1 for k“ and 11,.) the estimates of L' and I, were only slightly
improved: however. the estimate of 0,, was considerably more accurate. This would indicate that
the use of prior information had the greatest benefit on the estimation of (1,. especially when large
measurement enors were present. In all cases except the estimation of 11,, with measurement
enors with a large standard deviation and no prior information. the estimates of the thermal
properties were deemed acceptable. containing 1.3% error or less.
FromthesumofsquaresanalysisdescribedinCnapter4itwasexpectedthatthe
estimation of 0,, would require the fewest number of iterations for convergence. since its sum of
squares functionhadthe steepestminimum. Thesum ofsquares functions forL‘ and kdwere
much flatter. with the function for kd being slightly more steep than H. Therefore, it was
expected that the estimation of these properties would require more iterations. As shown in Tables
5.4, 5.5. and 5.6 for L'. k,,. and a, respectively the estimation of addid require the least mlmber
ofiterationswhenexactmeasurernemdatawasused,fouowedbyk,andL'. Theadditionof
randomenorsmthemeasurememdatadidnotalwaysresultinanimreasemmerequiredmunber
of iterations. The inclusion of prior information with a nail standard deviation did result in a
decreaseinthe requirednumberofiteraticns forallthreeproperties,especiallywhenthe standard
deviation of the measurement errors was large. As the standard deviation of the prior information
wasincreased.merequirednumberofiterafionshadatendencytoincrease. Thecaseoflalge
standard deviations for both the measurement errors and the prior information resulted in the

largest number of iterations required for convergence for all three thermal properties.

56

Table 5.4. Number of Iterations Required for Convergence in the Estimation of L‘

Standard Deviation of Measurement Errors

 

 

 

 

 

 

 

 

 

Table 5.5. Number of Iterations Required for Convergence in the Estimation of lg,

Standard Deviation of Measurement Errors
0.001 0.01 0.1

 

 

 

 

 

 

 

 

 

57
Table 5.6. Number of Iterations Required for Convergence in the Estimation of 01,,

Standard Deviation of Measurement Errors
0.001 0.01 0.1

 

 

 

 

 

 

 

 

 

5.1.2 Simultaneous Estimation of the Thermal Properties

Without the use of prior information the thermal properties L' and 01,, were simultaneously
estimated using exact temperature measurement data obtaimd from both t1: frozen and unfrozen
regions and with data containing random measurement errors. Prior information of the actual
valuesoftheseparameters wasthenincludedintheestimationproceduleandtheproperties were
again simultaneously estimated using exact data and data containing measurement enors. Three
different standard deviations for both the measurement errors and the prior information were used.
The results are presented in Table 5.7. Following the same procedure. the properties k. and 11,,

were also estimated simultaneously. with the results shown in Table 5.8.

58

Table 5.7. Simultaneous Estimation of the Dimensionless Latent Heat of Fusion. L‘
and the Dimensionless Thermal Diffusivity, a,

Standard Deviation of Measurement Errors

 

0.001

L’ 8 -99.998
1 0.013
%errcr 8 0.1112

0.01

L’ 8 -99.855
1 0.149
%error = 0.145

0.1

L. 8 -99.727
1 1.014
%error 8 0.273

U 8 400.001

%error 8 0.(I)1

 

a, 8 1.00109
1 0.00144
%error = 0.1119

n, 8 1.00053
1 0.00979
%arcr = 0.053

a, 8 1.09572
1 0.12598
%aror .-.- 9.572

11,8 1.001110

%errcr=0.m0

 

L. 8 -99.999
1 0.010
%a'ror 8 0.1111

L' 8 -99.989
1 0.011
%error 8 0.011

U 8 -99.999
1 0.000
%error 8 0.001

I: 8 -100.000

%error=0.(DO

 

a. = 0.99952
1 0.00115
‘56!!! 8 0.048

ad 3 1.111110
1 0.00003
%aror = 0.1110

(1,, = LIXXIX)
1 0.00000
%error = 0010

0,: 1.00000

%error=0.m0

 

L' 8 -99.999
1 0.013
%error 8 0.1111

L. 8 -99.872
1 0.130
%error 8 0.128

U 8 -99.992
1 0.082
%error 8 0018

L. 8 400.001

W = 0.001

 

a, 8 ”XXX”
1 0.00139
%crror = 0019

a, 8 1.00(Xl8
1 0.00230
%error = 0018

a, 8 1.00023
1 0.00039
%errcr = 0.023

(1,8 1.00000

meoooo

 

L. 8 -99.999
1 0.014
‘ %error 8 0.001

U 8 -99.855
1 0.150

' %error .. 0.14s

L. 8 39.878
:1: 0.963
m 8 0.122

U 8 -99.998

%err0r8 0.1112

 

 

a, 8 1.00009
1 0.00144
%error = 01119

 

a, 8 l.(XX)49
1 0.00948
%error 8 0.049

 

tr,' 8 1.01652
1 0.02891
%errcr 8 1.652

 

a, 8 ”XXX”
%error 8 0.001

 

59

Table 5.8. Simultamous Estimation of the Dimensionless Thermal Conductivity. 11,,
and the Dimensionless Thermal Diffusivity. (1,,

Standard Deviation ofMeasurement Errors

 

0.1!)!

k, 8 LW
1 0.00015
%error 8 0.1116

0.01

k, 8 0.99953
1 0.00173
%error 8 0.047

0.1

k, 8 1.01114
1 0.01459
%error 8 1.114

 

(1,, 3 0.99988
1 0.(Xl101
%error 8 0.012

or, = 0.99365
1 0.00859
%errcr = 0.635

a, 8 1.02600
1 0.08109
%error 8 2.600

 

kl a LW
1 0.00013
%aror 8 0.1114

k, 8 0.99998
1 0.00010
%error 8 0012

k‘ 3 Lam1
1 0.00010
%error 8 0.1111

 

11,8 0.99997

1 0.00024
%error 8 0.(X)3

a. 8 0.99998
1 0.00003
%error 8 0012

11,8 1.001110
1 0.001110
%error 8 0.1110

 

k‘ 8 ”XXIX
1 0.00015
%error 8 0014

k, 8 0.99970
1 0.00144
%errcr 8 0.030

k, 8 1.00066
1 0.00090
%error 8 0.1117

 

a, 8 0.99988
1 0.00098
%error 8 0.012

a, 8 0.99847
1 0.00204
%uror 8 0.153

(1,, 8 LW
1 0.00024
%error 8 0014

 

k, 8 ”KIDS
:1: 0.1XX115
m 8 0.1X15

k, 8 0.99953
1 0.00172
%error 8 0.047

k, 8 101926
:1: 0.01233
m 8 0.926

 

 

a, 8 0.99988
1 0.00101
%error 8 0.012

 

a. 8 0.99384
1 0.00832
%error 8 0.616

 

a, 8 1.00341
1 0.01863
%error 8 0.341

 

 

60
As shown in Table 5.7. the simultaneous estimation of L‘ and 01,, using exact measurement

data provided estimates that were highly accurate. containing only 0.001% and 0.000% error for
L‘ and a... respectively. The addition of measurement errors with standard deviations of 0.001,
0.01, and 0.1 to the estimation procedure resulted ill estimates that became less accurate as the
standard deviation increased. The presence oflarge measurement errors (a 8 0.1) provided an
estimate of L‘ containing only 0.273% error. but the estimate of (1,, was highly inaccurate.
comaining 9.572% error. The use of prior information with standard deviations of 0.1 and 0.001
for L' and 11,, respectively provided an overall increase in the accuracy of both estimates. In this
case. the simultaneous estimates obtained when large measurement enors (o 8 0.1) were present
were comiderably more accurate than those obtained without the use of prior information.
containing only 0.001% and 0.000% error for L’ and (1,, respectively. As the standard deviation
of the prior information was increased. the inaccuracy of the simultaneous estimates and their
corresponding 95% confidence linervals also increased. Using prior information with standard
deviations of 10.0 and 0.1 for L’ and 01,. respectively provided results that were more accurate than
those obtained without the use of prior information only when large measurement errors were
present. In this case, the estimates were considerably more accurate. containing 0.122% and
1.652% enor for L' and 01,, respectively.

The accuracy inthe simultaneous estimates oflcdand ordfollowed the same trend as
described above. as demonstrated in Table 5.8. Using exact measurement data in the estimation
procedure provided highly accurate estimates. containing only 0.001% and 0.000% error for kdand
a, respectively. The addition of random measurement errors with standard deviations of 0.001.
0.01,and0.1providedanoverallincleaseinboththeinacculacyoftheesfimatesandthe
corresponding 95% confidence intervals. Wlth large measurement enors (o 8 0.1) present. the
estimates of k“ and adcontain 1.114% and 2% error respectively. The addition of prior

information with a standard deviation of 0.001 for both parameters in the estimation procedure

61
resultedinanoveralldecreaseintheamountoferrorpresentinthesimultaneousestimates. As

the standard deviation of the prior information was increased to 0.01 and 0.1 for both parameters.
the inaccuracy ofthe estimates increased. The use ofprior information with a standard deviation
of0.1providedesfimatesthatwemmbederthandroseobtainedwidnfidreuseofpfior
information. exceptwhenlargemeasurementerrorswerepresent. Inthiscase,theestimatesof
1,, and 0,, were improved. containing 0.926% and 0.341% enor respectively.
lnthesimultaneousestimationofL‘andctp.itwasdeterminedthattheuseofprior
information with small standard deviations resulted in simultaneous estimates of significantly
higher accuracy than those obtained without prior information. especially when large measurement
enors are present. However, the use of prior information with large standard deviations had little
effect on the accuracy of the estimates, except when large measurement errors are present. In this
case. the estimate of L' was only slightly improved. but the estimate of a, was significantly
improved. When k,,and 12,, were simultaneously estimated. similar findings were obtained. Once
again, it was found that the use of prior information. even with a large standard deviation. had the
greatest benefitonthe accuracy ofthe estimateofct“ when large measurement enors were present.
With the exception of the estimate obtained for 01,, when the standard deviation of the
measurement enors was large and no prior information was used. all simultaneous estimates were
deemed acceptable with a maximum enor of less than 1.7%.
mbodlcases.dlemmberofitemfionsmquimdforconvergenceincreasedwimthe
addition of random measurement enors. as shown in Tables 5.9 and 5.10. The use of prior
information with small standard deviations resulted in an overall decrease in the required number
of iterations. especially when large measurement enors were present. As the standard deviation
of the prior information increased, so too did the number of iterations required for convergence.
The maximum required number of iterations occurred when both the prior information rmd

measurement error standard deviations were large.

62

Table 5.9. Number of Iterations Required for Convergence in the Simultaneous
Estimation of L’ and a,

Strmdard Deviation of Measurement Errors
0.001 0.01 0.1

 

 

 

 

 

 

 

 

 

Table 5.10. Number of Iterations Required for Convergence in the Simultaneous
Estimation of It, and 01,,

Standard Deviation of Measurement Enors
0.001 0.01 0.1

 

 

 

 

 

 

 

 

 

63
5.1.3 Comparison of Individual and Simultaneous Estimates

Both the individual and simultaneous estimates of the thermal properties. obtained with
the use of the Box-Kanemasu Interpolation Method, were extremely accurate with the exception
ofdnesfimwounmdforauusingtanperammmeasuremaudamwnmimnghrge
measurement errors (a 8 0.1) and no prior information. In this case the estimates were highly
inaccurate. containing as much as 9.6% enor. In all other cases the estimates were highly
accurate. containing less than 1.3% when estimated individually and 1.7% when estimated

simultaneously.

5.2 The Optimal Cryosurgical Treatment Time

As previously discussed. the accurate determination of the optimal cryosurgical treannent
time is of extreme importance to ensure adequate destruction of the cancerous tumor of a given
radius. while maintaining as much of the surrounding healthy tissue as possible.

5.2.1 Determination of the Optimal Cryosurgical Treatment Time using Prior Information
from the Sanre Radius

Without the use of prior information the optimal treatment time. 1,, was determined using
temperatureandexactradiusmeasurementdataobtaimdfrom boththefrozenandunfrozen
regions and with radius measurement data containing random measurement errors. Prior
information of the actual value of re. obtained from a previously performed procedure with the
sametumorradius, wasthenincludedintlreestimationprocedure;flletreatrnenttimewas again
dcterminedusingexaamdiusmeanlremaudamanddataconmimngrandomenors. Three
different standard deviations for both the radius measurement enors and the prior information were

used. The results are presented in Table 5.11.

64

Table 5.11. Estimation of the Optimal Treatment Time. 1.. using Simulated
Data from both the Frozen and Unfrozen Regions

Standard Deviation of Memuretnent Enors
0.0001 0.001 0.01

 

I, 8 0.18394 1‘ 8 0.17115
1 0.00074 1 0.00760
%error 8 0.573 %error 8 7.486

 

r. - 0.18416 1, - 0.17373 1. - 0.18497

1 0.00058 1 0.00601
%error 8 0.454 %error 8 6.092 %error 8 0.016

t. 8 0.18406 1‘ 8 0.17259 t‘ 8 0.18496
1 0.00065 1 0.00662
%error80.508 %errcr86.708 %aror80.022

r, 8 0.18406 1, 8 0.17258 1, 8 0.18496
1 0.00065 1 0.00663
%error 8 0.508 %ar0r 8 6.714 %error 8 0.022

 

 

 

 

 

 

 

The determination of 10 without prior information and with exact radius measurement data
resulted in an accurate estimated value. containing 0.022% error. as shown in Table 5.11. The
addition of random measurement enors with standard deviations of 0.0001. 0.001. and 0.01 to the
radius valuesresultedinanoveralldecreaseintheaccuracyoftheestimates,withanincreasein
the corresponding 95% confidence intervals. The estimate of 1, obtained using large measurement
errors (a 8 0.01) was highly inaccurate. comaining 7.486% error. The addition of prior
information with a standard deviation of 0.0001 provided only a slight increase in the accuracy
of the determined values of 1,. As the standard deviation of the prior information was increased
to 0.1111 and 0.01. there was very little change in the accuracy of the estimates. The value of 1,.
obtainedbyusingpriorinformationandradiusmeasurementerrors withstandalddeviationsof
0.01. contained 6.714% error. only slightly better than the value obtained without the use of prior

information. Except when large radius measurement enors were present. the determined values

65
for I, were deemed acceptable. containing less than 0.6% enor.

From the sensitivity coefficient study conducted prior to the estimation procedure. it was
determined that the most available information about the actual value of r, is contained in the data
obtained from the frozen region. Therefore, the optimal treatment time was again determined by
using temperature and radius meamrement data obtained from the frozen region only. with results

presented in Table 5.12.

Table 5.12. Estimation of the Optimal Treatment Time. tc. using Simulated
Data from the Frozen Region Only

       
 
 

Stmdard Deviation of Measurement Enors
0.0001 0.001 0.01

  

 

     

rt 8 0.18473
1 0.00004
%errcr 8 0.146

r, 8 0.18478

1, 8 0.18419
1 0.00063
%errcr 8 0.438

r, 8 0.18447
1 0.001113 1 0.00054 1 0.00428
%error80.ll9 %error80.236 %error8 3.865 %err0r80.114

r, 8 0.18476 tc 8 0.18442 1, 8 0.17716 tc 8 0.18477
1 0.001113 1 0.00059 1 0.00471
%error 8 0.130 %error 8 0.314 %error 8 4.238

t, 8 0.18476 re 8 0.18442 1, 8 0.17715
1 0.001113 1 0.00059 1 0.00471
W 8 0.130 W 8 0.314 %error 8 4.243

r, 8 0.17571
1 0.00554
%error 8 5.022

r. 8 0.17785

r. . 0.18477

    

 
 

 
 

    
     
  

%errcr 8 0.124
r. 8 0.18479

 

 

     

   

 
 

 
 

    

 

 

      
    

%error 8 0.124
r, 8 0.18477

      
 
 

   

 

    

  
 

 
 

    
   

%error8 0.124

 

 

 

 

  

Thistableshowsmatmeacalracyofmeesdmatesofrcexhibitsdlesamenendasthose

obtainedusinga-ulated measurementdatafromboththefrozenandrmfrozenregions. Itwas
expected fromthesensitivity coefficientanalysisthattheestimatesobtainedusingdatafromthe
frozen region only would be more accurate. However. when exact data was used without prior

information. the determined value of t, was less accurate. containing 0.124% enor. Also. the

66
estimates obtained using measurement errors with a standard deviation of 0.0001 and prior

mformafionwerelessaccumtethandloseobtahledusingdatafiombodrregions Asthestandard
deviation of the measurement errors was increased. the resulting estimates were more accurate than
those obtained using data from both regions. Except when large measurement errors (a 8 0.01)
were present. the estimates contain less than 0.4% error. With the presence of large measurement
errorstheestimatescontainasmuchas5%error.

It was also determined from the sensitivity coefficient analysis that very little information
aboutthe value oft. is contained intheunfrozen regiondata. However. theestimationprocedure
was repeated using temperature and radius measurement data obtained entirely within the unfrozen

region to determine if I, could be accurately estimated. The results are presented in Table 5.13.

Table 5.13. Estimation of the Optimal Treatment Time. 1,. using Simulated
Data from the Unfrozen Region Only

  
    

Standard Deviation of Measurement Enors
0.0001 0.001 0.01

  

 

     

        
  

I. 8 0.18432
1 0.00003
%error 8 0.368

t. 8 0.18310
1 0.00064 1 0.00619
%error 8 1.027 %error 8 8.438 %error 8 0.351

I, 8 0.18460 1‘ 8 0.18385 1, 8 0.17405 1‘ 8 0.18466
1 0.00004 1 0.00036 1 0.00344
%error 8 0.216 %error 8 0.622 %error 8 5.919 %error 8 0.184

t‘ 8 0.18435 1. 8 0.18321 1‘ 8 0.16994 re 8 0.18436
1 0.00032 1 0.00059 1 0.00471
%error 8 0.351 %errcr 8 0.968 %error 8 8.141 %error 8 0.346

r, 8 0.18434 1. 8 0.18320 1, 8 0.16988 1‘ 8 0.18435
1 0.00002 1 0.00059 1 0.00473
%error 8 0.357 %error 8 0.973 %aror 8 8.173

I, 8 0.16939 1‘ 8 0.18435

     
  

 
 

 
       
      

     

 

         
        

        
  

 

  

      
     
  

    
   

 

 

 

  
   

    

 
 

 
   

        
     

%error 8 0.351

 

 

 

 

  

 

As demonstrated in this table. the estimates of 1. exhibit the same tendencies as those

obtained using information from the frozen region only. However. there was an overall decrease

67

in the accuracy of the estimated values, as predicted by the sensitivity coefficient analysis. The
estimates contain less than 1.03% error except when large measurement errors (a 8 0.01) were
present. In this case. the estimates contain as much as 8.5% error.

Inmedetermmadonofdleopfimalneannaudmemqmredmachieveadesiredmmimum
temperatureataspecified location. itisconcluded thatthe use ofpriorinformation from the same
radius. regardless of the standard deviation. only slightly improved the accuracy of the estimates
of 1.. As expected from the sensitivity coefficient analysis. the estimates of r. obtained by using
data entirely within the frozen region provided the highest degree of accuracy. except when exact
dataordatacontainingmeasurernentenors withasrnall standarddeviationwereused. Inallcases
except when the large measurement errors (a 8 0.01) were present. the estimation prowdure
provided acceptable results. with estimated values of 1, containing 1.03% or less.

mememofsquamsmalysispresenmdmfllapter4.itwasexpeaeddmmenumber
of iterations required for convergence would be small due to the extremely steep minimum of the
sum of squares flulction. As demonstrated in Tables 5.14, 5.15, and 5.16. the number of iterations

remained mall and fairly constant.

Table 5.14. Number of Iterations Required for Convergence in the Estimation of 1.
using Simulated Data from both the Frozen and Unfrozen Regions

Standard Deviation of Measmement Errors
0.001 0.01

 

 

 

 

 

 

 

 

68

Table 5.15. Number of Iterations Required for Convergence in the Estimation of t,
usingSimulatedDatafromtheFlozenRegionOnly

Standard Dev'ution of Measurement Enors
0.0001 0.1111 0.01

 

 

 

 

 

 

 

 

Table 5.16. Number of Iterations Required for Convergence in the Estimation of 1,
using Simulated Data from the Unfrozen Region Only

Standard Deviation of Measurement Errors
0.0001 0.001 0.01

 

 

 

 

 

 

 

 

 

5.2.2 Determination of the Optimal Cryosurgical Treatment Time Using Prior Information
from a Different Radius

A more feasible source of prior information is obtained from a previously performed

procedure with a different tumor radius: therefore. the determination of the optimal treatment time

69
was repeated using prior information obtained from a single different radius location.
Measurement error and prior information standard deviations remained at 0.0001, 0.001. and 0.01.

The results are presented in Table 5.17.

Table 5.17. Estimation of the Optimal Treatment Time. 1,. using Prior
Information Obtained from a Different Radius

Standard Deviation of Measurement Enors
0.0001 0.001 0.01

1c 8 0.18394
1 0.00074
%error 8 0.573

1c 8 0.17115
1 0.00760
%mor 8 7.486

t, 8 0.17590
:1: 0.00745

 

t‘ 8 0.18404

1 011153

%error 8 0.519

%error 8 4.919

 

re 8 0.18402
1 0.00065
%error 8 0.530

r‘ 8 0.17579
1 0.00721
%error 8 4.978

 

 

 

I, 8 0.18401
1 0.00065
%error 8 0.535

 

tc 8 0.17577
1 0.00720
%error 8 4.989

 

 

As showninthistable.theuseofpriorinformationfrom adifferentradiuslocationwith
a standard deviation of 0.0001 provided a slight increase ill the accuracy of the estimates for
measurement enor standard deviations of 0.0001 and 0.001 when compared to the values obtained
without prior information. The estimate obtained using large measurement errors (a 8 0.01) was
considerably improved. containing 4.919% error. compared to an enor of 7.486% without the use
of prior information. An increase in the standard deviation of the prior information only slightly
decreasedtheaccuracyoftheestimates. Withthepresenceoflargemeasurementenomthe
estimates were again found to be highly inaccurate. containing as much as 5% enor.

To determine iftheaccuracyoftheestimates ofrccouldbeimproved.theestimation

procedure was repeated using prior information obtained from two previously performed

70
procedures with two different tumor radius locations. The results are shown in Table 5.18.

Table 5.18. Estimation of the Optimal Treatment Time, I... using Prior
Information Obtained from Two Different Radii

  
     

  

Standard Deviation of Measmement Errors
0.0001 0.001 0.01

 

     

     
      
    

      
  

t, 8 0.18486
1 0.00018
%error 8 0.076

r. 8 0.18483
1 0.00107
m 8 0.092

I, 8 0.18497
1 0.01118
%error 8 0.016

r, 8 0.18486
1 0.00017
%error 8 0.076

t‘ 8 0.18394
1 0.00074
%error 8 0.573

r, 8 0.18458
1 0.00050
%error 8 0.227

t‘ 8 0.18436
1 0.00057
%aror 8 0.346

r. 8 0.18432
1 0.00056
%error 8 0.368

1‘ 8 0.17115
1 0.111760
%error 8 7.486

I, 8 0.17279
1 0.00962
%error 8 6.600

t, 8 0.17243
1 0.00877
%error 8 6.795

t‘ 8 0.17123
1 0.01052
%error 8 7.443

r. 8 0.18496

 
 

   
 

   
 

%aror 8 0.022
r. 8 0.18495

 

      
 

 
 

 
  

 
  

%errcr 8 0.1D7
r. 8 0.18496

   
 
 

 

      
 

 
  

     
  

       

m 8 0.022
1‘ 8 0.18496

 

    

    
 

 
 

 
 

    
  

%errcr80.022

 

 

 

 

 

As demonstrated in this table. it was found that a small standard deviation (0 8 0.0001)

ofthepriorinformationresultedinanincreaseintheaccuracyoftheestimatesoft,.exceptwhen
thestandarddeviationofthemeasurementerrorswassrnall. Inthiscase,theuseofprior
information from two different radii actually resulted in a slight decrease in the accuracy of the
estimated value of 1.. As the standard deviation ofthe prior information increased. the accuracy
oftheestimatesdecreasedslightly. Withitllepresenceoflargemeasurementerrors(o80.01),
theestimateswerefoundtocomainasmuchas75%enor. Inallothercases.theestimateswere
formd to be highly accurate. containing less than 0.4% enor.
Indledetermmadonofmeopdmalueamentdmeusingpmrmfomafionobtainedfiom
a single different radius location, the number of iterations required for convergence was again
foundtobesrnanandfainyconstantasshowninTable519. Thenumberofiterationsrequired

for convergence when prior information from two different radius locations was used was also

71
found to be small and fairly constant. with the exception ofthe case of large standard deviations

of both the radius measurement errors and the prior information. These results are shown in Table

5.20.

Table 5.19. Number of Iteratiom Required for Convergence in the Estimation of I.
using Prior Information Obtained from a Different Radius

Standard Dev'ntion of Measmement Enors

 

 

 

 

 

 

 

 

Table 5.20. Number of Iterations Required for Convergence in the Estimation of 1,.
using Prior Information Obtained from Two Difierent Radii

    

Standard Deviation of Measurement Enors
0.11101 0.001 0.01

  

 

     

 

 

 

 

 

 

 

72

From this portion of the investigation. it was found that the use of prior information
obtained from a previously performed procedure with a single difierent radius provided estimates
of t. that are only slightly more accurate than those obtained without prior information. The
exception to this is when large measurement errors were present in the estimation prowdule. 111
this case. the use of prior information from a different radius location provided considerable
improvementintheaccuracyoftheestimate. 'Iheuseofpriorinformationfromtwodifi’erent
radiuslocationsdidnotprovideasignificantincreaseintheaccuracyoftheestimates. Itis
comludedthatdreuseofpriorinformationobtainedfiomasingleradius locationisquite

beneficial when large radius measurement errors are present.

CHAPTER 6

SUMMARY AND CONCLUSIONS

The primary goal of this invenigation of the cryosurgical freezing of undesirable tissue
was to test the minimization procedure. the Box-Kanemasu Interpolation Method. for accuracy and
reliance in the estimation of the tissue thermal properties and the determination of the optimal
treatment time required to achieve a desired minimum temperature at a specified radius location.
The estimated values provided by this procedure were compared with the actual values used to
generate the simulated meanlrement data required as input for the Box-Kanemasu Method. It was
found that the methodologies presented in this study provided very accurate estimates of the

thermal properties and optimal cryosurgical treatment time.

6.1 Estimation of the Thermal Properties

The estimation of the thermal properties. namely the dimensionless latent heat of fusion.
L‘, the dimensionless thermal conductivity. It, and the dimensionless thermal diffusivity, a... was
wnductedusmgexaatempemnuemeamrunandaamddmawmainmgmndanmeasummem
errors. both with and without prior information. The results of this portion of the investigation
support the following conclusions:
1. Both with and without the use of prior information the Box-Kanemasu Interpolation Method

provided extremely accurate estimates of the thermal properties with errors of less than

73

74
1.3% when estimated individually and 1.7% when estimated simultaneously. The only

excepdonwasdmesfimatesobtainedfora,whenmeasuremahenom%astmdard
deviation. 0. of 0.1 (~ 1%) were present. In this case. errors as large as 9.6% were prwent.
The use of prior information provided an overall increase in the accuracy ofthe estimated
thermalproperties. Pdorinformafionhaddregreatesteffectontheesfimafionofmwhen
large measurement errors were present by considerably reducing the amount of error in the
estimated values. Therefore. prior information. if available. should be included in the

estimation procedure. particularly in the estimation of ad with large measurement errors.

6.2 Determination of the Optimal Cryosurgical Treatment Time

The determination of the optimal treatment time was performed using exact radius

measurement data and with data containing random measurement enors. both with and without

prior information obtained from a previously performed procedure with the same tumor radius.

The estimation procedure was repeated using prior information obtained from a previously

performed procedure with both one and two different tumor radii. The following conclusions are

supported by this portion of the investigation:

1.

Without the use of prior information. the Box-Kanemasu Interpolation Method provided
highly accurate estimates of the optimal treatment time, except when large measurement
enors (o 8 0.01. or ~10% ) were present In this case, enors as high as 7.5% were present.
The use of prior information obtained from the same radius location offered little
improvement in the accuracy of the estimated values of 1., for all standard deviations of
meamrement enors.

When large measurement enors were present. the use of prior information obtained from

anngledifieranmdiuspmvidedsignificnumcmasemtheaccuracyofdwesfimated

75
values. while the use of prior information from two radius locations resulted in only slight

improvement. Therefore. the use of prior information from a single different radius location
shouldbeusedinthedeterminationoftc. especially withthepresenceoflarge radius

11168811131110!!! CI'I'OI'S.

CHAPTER 7

RECOMMENDATIONS FOR FUTURE WORK

The methodologies presented in this investigation were found to provide very accurate
estimates of the thermal properties and optimal cryosurgical treatment time. However. because
the problem had been simplified. such as assuming the tissue to be homogeneous with constant
thermal properties. the actual freezing process of the tissue subject to cryosurgical freezing is not
accurately described. It is suggested that non-constant thermal properties be incorporated into the
procedure. These properties would be both temperature and spatially dependent to account for
thetruenanlreofthetissue. Also.irregularshapedgeometriesneedtobeconsidered.since
malignant tissue is rarely spherical in shape. These changes would require the use of a numerical
method to solve the describing differential equations. such as a finite element analysis.

It is also felt that the heat transfer equations do not adequately describe the heat transfer
process of the tissue before and during freezing. From the Literature Review presented in Chapter
2.itwasfoundthatdreefl’ectsofmebloodperfusionmteonthefieezingpmcesswere quite
significant. Also. the effects of the metabolic heat generation rate, while not as great. would need
to be investigated further. Therefore. it is recommended that the bieheat equation. equation (2.1).
beusedtodescribetheheatuansferirlthelmfiozenregionofthetissue,withtheheatequation
used for the frozen region.

Experimental work is also required to further test the estimation procedures for accuracy.

Thisworkcouldbeconductedonsimulatedtissue.suchasthegelatinsolutionpresentedin

76

77
Chapter 2. or on laboratory animals. In both cases. the temperature history of the tissue being

frozen would be recorded and used as input for the Box-Kanemasu Interpolation Method in lieu
of the simulated measurement data. Another source of measurement data could be obtained from
a previously performed procedure in which thermocouples were used as the monitoring device.

This procedure could be extended to include the determination of the optimal location of
thecryoprobe withinthetissuetobeflozen. Asthetumorisnolongerconsideredtobe
homogeneous and is considered to be irregular in shape. the location of the probe will not be at
the geometric center. The radius of the tumor would be given as

r -- [(x-x’)’ + 0-y’)’ + (z—z’fl‘” (7.1)

In this case. the objective flmction would be minimized with respect to the probe location.
(x’.y’.z').

In this investigation. the presented methodologies have proven to be accurate and reliable
in the estimation of the thermal properties and the determination of the optimal cryosurgical
treatment time. Therefore. this work is a solid foundation upon which to build and expand the

aforementioned recommendations.

APPENDICES

APPENDIX A

THE FORTRAN PROGRAM SENSE.F OR

This program. SEN SE.FOR. is used to calculate the sensitivity coefficients for the thermal

properties to be estimated.

PROGRAM SENSE

THIS PROGRAM IS DESIGNED TO CALCULATE THE SENSITIVITY

COEFFICIENTS OF THE TEMPERATURES,WITH RESPECT TO L.

KSL, AND ALPHASL. FOR BOTH THE FROZEN AND UNFROZEN

REGIONS SURROUNDING A POINT SOURCE HEAT SINK.
WRITTEN BY LESLIE SCO'IT

0000000

DOUBLE PRECISION KSL, Q. ALPHSL. L. PI
DOUBLE PRECISION ERFC. ZBRENT
DOUBLE PRECISION DETA. ALPHSR. PISR. X. XX. X2. X2ALPH.

+ XALPH. XXALPH, ET A, ETAS. ETA2. BO’I'I‘OM, TOP.
4» TOP2, DER]. DER2A, DER2, DER3, DER4, DERS,
+ DER6. DER7, DER8. SENSE]. SENSEZ. SENSE3.
+ SENSFA, DSENSEI, DSENSEZ. DSENSE3. DSENSE4
C
COMMON/PROP/KSL. Q. ALPHSL. L. PI
C
EXTERNAL ERFC. ZBRENT
C
OPEN(UNIT=10. FILE="SENSE.DAT". STATUS="UNKNOWN”)
L '5 -1(X).D0
Q = -l.DO
KSL = l.D0
ALPHSL = 1.D0
P1 = DACOS(-1.DO)
NT = 2(1) '
INCR = l

DETA = 1.0D—2

78

0 00000

000000000000000000000000

VARIABLES DECLARED

ALPHSR = A1.PHSL*"'O.5DO
PISR = (PI“0.SDO)f2.DO

X IS THE CALCULATED VALUE FOR LAMBDA
ZBRENTisamotflndingsubmufineusedmsolvetheuanscmdentalequafionforme
freezing front location See Numerical Recipes by Press et al.. Cambridge University
Press. New York. New York. 1986.

X = ZBRENT(1.0D4.2.D0.0.00IDO)
WRIT‘E(10.*)”X=",X
XX = X‘X
X2 = 2.D0*X
X2ALPH = X2*ALPI-ISR
XALPII = X‘AIPHSR
XXALPH = XX'ALPHSL
DO I = 2.NTJNCR
ETA = (I-l)‘DET‘A
ETAS = BTA*ETA
ETA2 = 2.0DO*ETA
BOTTOM = (EXP(-XXALPH))/X2ALPH - PISR'ERFC(XALPH)
TOP = (EXP(-XXALPH))/(2.0D0*XX*ALPHSR)
TOP2 = (EXH-AIPHSL’ETASMET‘ArALPIISR» - PISR
+ *ERFC(ETA*ALPHSR)

DER] = THE DERIVATIVE OF THE TEMPERATURE FUNCTION. IN THE
SOLID REGION. WITH RESPECT TO LAMBDA

DERZA = THE DERIVATIVE OF TIE FUNCTION USED TO SOLVE FOR

LAMBDA WITH RESPECT TO LAMBDA
DER2 = TIE INVERSE OF DERZA
DER3 = THE DERIVATIVE OF TIE FUNCTION USED TO SOLVE FOR
LAMBDA WITH REPECT TO L

DER4 = TIE DERIVATIVE OF TIE FUNCTION USED TO SOLVE FOR
LAMBDA WITH RESPECT TO Q. NOT USED

DERS = TIE DERIVATIVE OF TIE FUNCTION USED TO SOLVE FOR
LAMBDA WITH RESPECT TO KSL

DER6 = TIE DERIVATIVE OF TIE FUNCTION USED TO SOLVE FOR
LAMBDA WITH RESPECT TO ALPHASL

DER? = TIE DERIVATIVE OF TIE TEMPERATURE FUNCTION. IN TIE
LIQUID REGION. WITH RESPECT TO LAMBDA

DER8 = TIE DERIVATIVE OF TIE TEMPERATURE FUNCTION. IN TIE
LIQUID REGION. WITII RESPECT TO ALPHASL

SENSE] = TIE SENSITIVITY COEFFICIENT FOR TIE TEMPERATURE
WITH RESPECT TO L
SENSE2 = TIE SENSITIVITY COEFFICENT FOR TIE TEMPERATURE
WITH RESPECT TO Q. NOT USED
SENSE3 = TIE SENSITIVITY COEFFICENT FOR TIE TEMPERATURE

0000000

C
C

80

WITH RESPECT TO KSL
SENSE4 = TIE SENSITIVITY COEFFICENT FOR TIE TEMPERATURE
WITH RESPECT TO ALPHASL

CALCULATION OF TIE SENSITIVITY COEFFICENTS

DER] = Q‘(—EXP(-XX)/(2.0D0"XX))

DERZA = KSL*Q*((-XX*EXP(-XX) - EXH-XX))/(X*X‘X)) -

+ (((-XX*EXP(-XXALPH) - EXP(-XXALPI'I))/(X"‘X"'X"l

+ ALPHSR))*BOTT OM - TOP‘(((-2.0D0*XX"ALPHSL“

+ (3.0DOIZ.OD0)*EXP(-XXALPH) - ALPHSR‘EXP

+ (~XXALPH))/(2.0DO"XX*ALPHSL)) + ALPHSR*EXP

+ (-XXALPH))/(BOTTOM*BOTTOM)) - L

DER2 = -DER2A‘""(-].ODO)

DER3 = -X

DER4 = KSL‘EXPGXXMZDDO‘XX»

DER-5 = Q*(EXP(-XX)/(2.0DO*XX))

DER6 = o(((-2.0DO*XX*ALPHSR*EXP(-XXALPH) - ALPHSL

"'*(-0.5D0)“EXP(-XXALPH))/(4.0U)*XX*ALPHSL))

‘BOTT OM - TOP*((-2.0D0*XX*ALPHSR*

EXP(-XXALPH) - AI..PHSL“(-O.5D0)*EXP

(-XXALPH))/(4.0D0*X*ALPHSL) + 0.5D0*ALPHSL

""(-0.5D0)*X"EXP(-XXALPH)))/(BOTTOM*BOTTOM)

7 = -TOP2*((-2.0DO*XX*(ALPHSL**(3.0D0/2.0D0))*EXP

(-XXALPH)) - (ALPHSR‘EXP(-XXALPH))/(2.0DO*XX*ALPHSL)

+ ALPHSR*EXP(-XXALPH))/(BOTTOM*BOTTOM)

DER8 = ((((-2.0D0*(ETA“3.0D0)"ALPHSR*EXP(-ETAS*ALPHSL))

-(EXP(-ETAS"ALPHSL)‘ETA*(ALPIISL“(—O.5D0))))

l((2.0DO*ETA*ALPHSR)“2.0DO) + (0.5DO*ALPHSL

“(-O.5DO)"'ET A*EXP(-ETAS*ALPHSL)))“BOTT OM

- P2*((-EXP(-XXALPH))*(XX + (2.0D0‘X‘ALPHSL

*DER2‘DER6))*(2.0D0*XALPH) - ((EXP(-XXALPH))*
((2.0D0‘AIPHSR‘DER2‘DER6) + (X‘ALPIISL“(-O.5D0))))/

((2.0DO‘XALPH)"2.0D0) + (EXP(-XXALPH)*(O.5D0“‘
ALPHSL“(-0.5D0)"X + AIPHSR‘DERZ‘DER6))))

I(BOTTOM*BOTTOM)

++E+++++

+++++++++

CALCULATION OF FROZEN PORTION SEN SITIVTTY COEFFICENTS
IF(ETA .LT. X) TIEN
SENSE] = DER]‘DER2*DER3
DSENSE1= L‘SENSE]
SENSE2 = -(EXP(-ETAS)/ETA2 -- EXP(-XX)/X2 - PISR

+ *(ERFCGETA) - ERFCOO» - Q“(DER1"
+ DER2‘DER4)
DSENSE2=Q*SENSE2
SENSE3 = DER1*DER2*DER5
DSENSE3=KSL*SENSE3

SENSE4 = DERI‘DERZ‘DER6

8]

DSENSE4=ALPHSL*SENSE4
ELSE
CALCULATION OF UNFROZEN PORTION SENSITIVITY
COEFFICENTS
IF(ETA .GT. X) TIEN
SENSE] a DER7‘DER2‘DER3
DSENSE1=L"SENSE1
SENSEZ = DER7‘DER2‘DER4
DSENSEL-Q‘SEN SE2
SENSE3 = DER7‘DER2‘DER5
DSENSE3=KSL"SENSE3
SEN SE4 = DERS
DSENSE4=ALPIISL*SEN SE4
ELSE
SENSITIVITY COEFFICENTS AT TIE INTERFACE
SENSE] = 0.0D0
SEN SE2 = 0.0D0
SENSE3 = 0.0D0
SEN SE4 = 0.0D0
ENDIF
ENDIF
WRITE(10.5)ETA. DSENSEI
5 FORMAT(E]2.5.E]2.5)
ENDDO
STOP
END
CALCULATION OF LAMBDA FROM FUNCTION ZBRENT
DOUBLE PRECISION FUNCTI ON ZBRENT(X] . X2, TOL)

PARAMETER(ITMAX=lm. EPS= 3.]E-8)
DOUBLE PRECISION A, B, C. D. E. FA. FB. FC
DOUBLE PRECISION TOL], TOL. X1, X2. XM
DOUBLE PRECISION P. Q. R. S. FUNCL

EXTERNAL FUN CL
A=X1

B=X2
FA=FUN(1(A)
FB=FUNCL(B)

IF(FB"FA .GT. 0.0D0) PAUSE ’ROOT MUST BE BRACKETED FOR
+ ZBRENT.’
FC=FB
DO 15 ITPR=IJTMAX
IF(FB*FC .GT. 0.000) THEN
C=A
FC=FA
D=B-A
E=D

ENDIF

82

IF(ABS(FC) .LT. ABS(FB)) TIEN

=B
B=C
C=A
FA=FB
FB=FC
FC=FA
ENDIF

TOL]=2.0DO*EPS*ABS(B)+0.5DO‘TOL

XM=0.5DO‘(C-B)

IF(ABS(XM) .LE. TOL] .OR. FB .EQ. 0.0D0)TIEN

ZBRENT=B
RETURN
ENDIF

IF(ABS(E) .GE. TOL] .AND. ABS(FA) .GT. ABS(FB» TIEN

S=FBIFA

IF(A .EQ. 0mm
P=2.0DO*XM*S
Q=l.OD0 - s

ELSE
Q=FAIFC
R=FBIFC

P=S"'(2.0D0"XM*Q*(Q-R) - (B-A)*(R-l.OD0))
Q=(Q-1.0DO)*(R-I.ODO)*(S- 1.0D0)

ENDIF

IF(P .GT. 0) Q = -Q
P=ABS(P)

IF(2.0D0*P .LT. MIN(3.0D0'*XM*Q-ABS('TOL1‘Q).ABS(E*Q)))TIEN

E=D
D=P/Q
ELSE
D=XM
E=D
ENDIF
ELSE
D=XM
E=D
ENDIF
A=B
FA=FB
IF(ABS(D) .GT.T‘OL1)T'I-IEN
B=B+D
ELSE
B=B+SIGN(TOL1.XM)
ENDIF
FB=FUNCL(B)

83

15 CONTINUE

C

PAUSE ’ZBRENT EXCEEDING MAXIMUM TI'ERATIONS.’
ZBRENT=B

RETURN

END

DOUBLE PRECISION FUNCTION ERFC(X)

DOUBLE PRECISION A1. A2. A3. A4. A5. P. T, X
Al=0.254829592D0
A2=-0.284496736D0
A381.42]4]3741D0
A4=- l .45315202‘7D0
A5=].061405429D0
P=O.327591 1DO
T=l .ODO/(l .0D0+P"X)
ERFC=(A1*T+A2*T“2.0D0+A3"T"3.0DO+A4‘T“4.0DO+A5*T“5.0D0)
+ ‘EXP(-X“2.0DO)
RETURN
END

DOUBLE PRECISION FUNCTION FUNCL(X)

DOUBLE PRECISION KSL. Q. ALPHSL. L, X. P]. ERFC
DOUBLE PRECISION EXPX. EXPXASL. XX2. RATIO
COMMON/PROP/KSL. Q. ALPHSL. L. PI

EXTERNAL ERFC

EXPX=EXP(-X*X)

EXPXASL==EXP(-X*X*ALPHSL)

XX2=X*X*2.0D0
RATIO=EXPXASIJ(XX2‘ALPHSL“0.5D0)
FUNCL=KSL"Q*EXPX/XX2 - RATIO/(RATIO—(PI“0.5D0/2)
+ *ERFC(ALPHSL”O.5DO*X)) -L"'X

RETURN

END

84
This file represents the output file from the program SENSEFOR. The first column is

the dimensionless similarity variable I]. the second column contains the dimensionless sensitivity

coefficient values for L'.

.ICXXDEcO] -.44089E+m
.ZIXDOE-Ol -.44089E+(X)
.300005-01 -.44089E+(X)
.WE-O] -.44089E+(X)
WED] -.44089E+(X)
.soooora-or -.44089E+(X)
.7(Xl)OE-Ol -.44089E+(X)

.soooos-or ~.44089E+(X)
900005-01 -.44089E+(X)

. ](X)00E+00
.1 1000E+00
.]2(X)0E+00
. ] NEW
. “(HEAD
.ISIXXIE+(X)
. 160CDE+(X)
. I7000E+00
. 18000E+00
. 19ooor=.+oo
.20000E+00
.2 1(XX)E+(X)
flmOE-t-(X)
230305“)
.24OCDE+(X)
WEI-00
.ZGIXXEHX)
.27000E-t-(X)
.28(DOE+00
.29000E+00
.WEH'D
.31(XX)E+(X)
.32(IX)E+(X)
.33(IX)E+(X)
.3aooor=.+oo
.35000E-t-(X)
.mEt-(X)
.37ooor-:+oo
SW
.39IXX)E+(X)
.W
.4](XX)E+(D
AMER-(X)
.4aooor=.+oo

-.44089E+00
-.44089E+00
-.44089E+(X)
-.44089E+(X)
-.44089E+(X)
-.44089E+(X)
-.44089E+(X)
-.44089E+(X)
-.18875E+(X)
-. l7529E-t-(X)
-. 16322E+CD
-.]52345+00
-. 14250E-l-(X)
-.133SSE+(X)
-.]2538E+(X)
-.1 1790mm
-.1 1 104mm
-.10471E+00
-.98867E—01
-.93458E-0]
-.88439E-0]
-.83771E-0]
-.79422E-01
-.75363E-0]
-.71568E-01
-.68014E-01
-.64681E-O]
-.6]550E-0]
-.58606E-0]
-.5583SE-01
-.53222E-Ol
~50757E-0]
-.48428E-01
-.46226E-0]

.44(XX)E+(X)
AME-ta)
.46000E+00
.47(X)0E+(X)
.48(X)0E+(X)
.49000E-I-00
WEI-(X)
.5 ](X)0E+(X)
.SZIXIOE-t-(X)
53WE+IXI
.54000E+00
550009.00
560(1)E+(X)
.57000E+00
.58IXX)E+00
.59000E+00
.6(XX)0E+(X)
£10me
.62000E+(X)
.63m0E-t-(X)
.64000E+00
.6SOOOE+00
.66(X)0E+(X)
.67000E+(X)
.68000E+00
.69000E+00
.7(X)00E+00
.71(X)0E+(X)
.72(X)OE+(X)
.73(X)OE+(X)
.74000E+00
.75(X)0E+00
.76(X)0E+00
.77IX)OE+OO
.7sooor=.+oo
.79000E+00
.8(XXXIE+(X)
.8 100051.00
.82(XX)E+(X)
.83(XX)E+(X)
.84000E4-CX)
£50me
.86(X)OE+CD
.87(X)0E+(X)
.88(XX)E+(X)
.890(X)E+(X)
.9000054-00
.9rooor-:+oo
.9200054-00

-.44142E-Ol
-.42168E-Ol
-.40297E-Ol
-.38522E-01
-.36837E-01
-.35236E-01
-.337 l4E-Ol
-.32266E-Ol
-.30889E-Ol
~29576E-01
-.28326E-01
~27] 34E-O]
-.25998E-Ol
-.24913E-01
-.23878E-Ol
-.22889E-O]
-.21944E-01
-.2]042E-01
-.20]79E-Ol
-. I9354E-01
-. l 8564E-01
-. l7809E-0]
-. ] 7086E-01
-. l 6393E-Ol
-.lS730E-01
-. 150955-01
-.14487E-O]
-. ] 3904E-O]
-.] 3345E—01
-.12809E-0]
-.12295E-01
-.1 1802E-01
-.l l330E-01
-. 10876E-01
-. 10441E-Ol
-.l(XJ24E-Ol
-.96230E-02
-.92384E-02
-.88693E-02
-.85 1495-02
-.8 I746E-02
-.78479E-02
-.75342E-02
-.72329E-O2
-.69435E-02
-.66655E-O2
-.63986E-0Q
-.6142 lE-OQ
-.5895 8E-02

85

APPENDIX B

THE FORTRAN PROGRAM SQUARES.FOR

This program. SQUARESFOR. is used to calculate the sum of squares function for each

thermal property to be estimated.

PROGRAM SQUARES
C
C THIS PROGRAM IS WRITTEN TO CALCULATED THE SUM OF SQUARES
C FUNCTION USING EXACT DATA FOR Ya) AND DATA WITH A SINGLE
C PARAMETER VARIED FOR T(I)
C WRITTEN BY LESLIE A. SCOTT

DOUBLE PRECISION KSL. Q. ALPHSL. L. PI, DALPHSL
DOUBLE PRECISION DETA. BETA. TIE'TA
DIMENSION Y(5(X)). T(5(X))

COMMON/PROP/KSL. Q. ALPHSL. L. PI
COMMON TIETA. BETA. I

OPEN(UNIT=]O. FILE=”SQUAR.DAT". ST ATUS="UNKNOWN")

KSL a 1.000

Q = -l.OD0

ALPHSL = 1.000

L = -l(X).OD0

PI = DACOS(-].OIX))
NT = 150

INCR = I

DETA = 1.0D-2

DOI=2.NT.INCR

CALL MODEL
Y0) = TIETA

86

0

0 0000

87
ENDDO

ALPHSL = 0.50DO
DALPHSL = 0.050D0
N = 2]

DO I = U"
DO I = 2.NT.INCR
CALL MODEL
T0) = TIETA
ENDDO
SUM = 0.0D0
DO I = 2.NT.INCR
SUM = (Y (I) - T0))“(Y(I) - TO»
ENDDO
WRITE(]O.5)ALPIISL. SUM
5 FORMAT(EIZ.5.I312.5)
SUM = 0.000
ALPHSL = ALPHSL + DALPHSL
ENDDO
STOP
END

SUBROUTINE MODEL

DOUBLE PRECISION KSL, Q. ALPHSL. L. P]

DOUBLE PRECISION ERFC. ZBRENT

DOUBLE PRECISION DET A. ALPHSR, PISR. X. XX, X2. X2ALPH.
+ XALPH. XXALPH. BETA. BETAS. BETA2. TIETA

COMMON/PROP/KSL. Q. ALPHSL. L. PI
COMMON TIETA. BETA. I
EXTERNAL ERFC. ZBRENT

NT = 150
INCR = 1
DETA 81.0D-2

ALPHSR = ALPHSL“0.SDO

PISR a (PI“O.SDO)/2.D0
X 18 THE CALCULATED VALUE FOR LAMBDA
ZBRENTisamotfindingsubmudneusedwsolvemeuanscendemalequafionforme
freezing from location. See Numerical Recipes by Press et al.. Cambridge University
Press. New York. New York. 1986.

X = ZBRENT(1.0D4.2.D0.0.(X)1D0)
WRITE(10,*)"X=".X

XX = X‘X

X2 = 2.DO*X

X2ALPH = X2*ALPHSR

+

XALPH = X‘ALPHSR
XXALPH = XX‘ALPIISL
CALCULATION OF DIMENSIONLESS TEMPERATURES
BETA = (I-1)*DETA
BET AS = BET A‘BETA
BET A2 = 2.D0*BETA
CALCULATION OF FROZEN PORTION TEMPERATURE
IF(BETA .LT. X) TIEN
TIETA = 1-Q*(EXP(-BETAS)/BETA2 - EXP(-XX)/X2
-PISR*(ERFC(BETA) - ERFC(X)»
ELSE
CALCULATION OF UNFROZEN PORTION TEMPERATURE
IF(BETA .GT. X) TIEN
TIETA = (EXP(-ALPHSL"BETAS)/(BETA2*ALPHSR)
- PISR‘ERFC(ALPHSR*BETA))/(EXP(-XXALPH)/X2ALPH
-PISR*ERFC(XALPH))
ELSE
TEMPERATURE AT TIE INTERFACE. DETERMINED FROM B.C.
TIIETA = ].D0
ENDIF
ENDIF
WRITE(10.’(I10.7F10.5)’)I-1.TIETA.STDDV.BETA
RETURN
END
CALCULATION OF LAMBDA FROM FUNCTION ZBRENT
DOUBLE PRECISION FUNCTION ZBRENT (X1. X2. TOL)
PARAMETER(ITMAX=1(X). EPS= 3.0E-8)
DOUBLE PRECISION A, B. C. D. E, FA. FB. FC
DOUBLE PRECISION TOL], TOL. X1. X2. XM
DOUBLE PRECISION P, Q. R. S. FUNCL
EXTERNAL FUNCL
A=X1
B=X2
FA=FUNCL(A)
FB=FUNCL(B)
IF(FB*FA .GT. 0.0D0) PAUSE ’ROOT MUST BE BRACKETED FOR
ZBRENT.’
FC=FB
DO 15 TIER=].ITMAX
IF(FB*FC .GT. 0.0D0) TIEN
C=A
FC=FA
D=BoA
E=D
ENDIF
IF(ABS(FC) .LT. ABS(FB)) TIEN
=B
B=C
@A

89

FA=FB
FB=FC
FC'-FA
ENDIF
TOL1=2.0D0"EPS"ABS(B)+O.5DO"TOL
XM=0.5D0"(C-B)
IF(ABS(XM) .LE. TOL] .OR. FB .EQ. 0.0D0)TIEN
ZBRENT=B
RETURN
ENDIF

IF(ABS(E) .GE. TOL] .AND. ABS(FA) .GT. ABS(FB)) TIEN
S=FB/FA
IF(A .EQ. C)TIEN
P=2.0D0*XM*S
Q=1.0D0 - S
ELSE
Q=FA/FC
R=FBIFC
P=S*(2.0D0"XM‘Q*(Q-R) - (B-A)‘(R-I.OD0))
Q=(Q-].0DO)*(R-].OD0)‘(S-].OIX))
ENDIF
IF(P .GT. 0) Q = -Q
P=ABS(P)
IF(2.0D0"'P .LT. WN(3.0D0"XM*Q~ABS(TOL]*Q),ABS(E*Q)))TIEN
E=D
D=P/Q
ELSE
D=XM
ED
ENDIF
ELSE
D=XM
ED
ENDIF
ABE
FA=FB
IF(ABS(D) .GT. TOLI)THEN
=B+D
ESE
B=B+SIGN(TOL].XM)
ENDIF
FB=FUNCL(B)
CONTINUE

PAUSE ’ZBRENT EXCEEDING MAXIMUM I'TERA'IIONS.’
ZBRENT=B

RETURN

END

90

C
DOUBLE PRECISION FUNCTION ERFC(X)
C
DOUBLE PRECISION Al, A2. A3. A4. A5. P. T. X
Al=0.254829592D0
A2=-0.284496736D0
A3=1.42]4l374lD0
A4=-l.453]52m7D0
A5=1.06]405429D0
P=0.327591 IDO
T= I .0D0/( l .0DO+P"'X)
ERFC=(A1*T+A2*T"*2.0D0+A3"T“3.0D0+A4*T“4.0D0+A5“T“SDDO)
+ ‘EXP(-X“2.0DO)
RETURN
END
C
DOUBLE PRECISION FUNCTION FUNCL(X)
C
DOUBLE PRECISION KSL. Q. ALPHSL. L, X. PI. ERFC
DOUBLE PRECISION EXPX. EXPXASL. XX2. RATIO
COMMON/PROP/KSL. Q. ALPHSL. L, PI
EXTERNAL ERFC
EXPX=EXP(-X*X)
EXPXASL=EXP(-X*X*ALPHSL)
XX2=X*X*2.0DO
RATIO=EXPXASIJ(XX2*ALPHSL“O.5DO)
FUNCTfiKSL‘Q‘EXPX/XXZ - RATIO/(RATIO-(PI**O.5D0/Z)
+ *ERFC(ALPHSL“O.5DO*X)) -L*X
RETURN
END

APPENDIX C

THE FORTRAN PROGRAM NLINA.F OR

This program. NLINA.FOR. uses the Box-Kanemasu Interpolation Method to estimate the

thermal properties without prior information.

PROGRAM NLINA
CCCCCCCCC PROGRAM DESCRIPTION CCCCC
C C
C PROGRAM NLINC
C WRITTEN BY JAMES V. BECK . C

C LAST REVISED MAY 1. 1991
C REVISED BY LESLE A. SCOTT

C
Ct«utar-nustoutur-ueta-ususeteeter-«cutest"usueuuuuuuc
C C
CV CCCCCCCC VARIABLE IDENTIFICATION CCCCCCCC
C C
C C
centers-unu-«a:cueuteenintersects-eunu-sea-useetueuuuuuuc
C C
CDCCCCCCCC DIMENSION BLOCK BLOCK (XXX)
C C
C C
IMPLICTT REAL‘8 (A-H.O-Z)
DIMENSION T(35(X).5).Y(35W).SIGZ(35(X)).B(5).Z(5),A(5).BS(5).
IVINV(5.5).BSS(5).m(5),BSV(5).R(5.5).EXTRA(20),ERR(35(X))
1. PS(5.5).P(5.5).PSV(5.5).
I XT X(5.5).XTY (5)
CHARACTERMO DFILE.OUTFIL
C C
caucuses-nuueta-eueeecuecusect"cue“oneecuueuuuucuuc
C C
COCCCCCCCC COMMON BLOCK BLOCK 0100
C C

COMMON SIG2.TLBSJETAPSPBA.Y.MODL.VINV.NP.EXTRA

91

COMMON/ERROR/ERR
C C
C C
Cue«auauasueerr-eecueteeeeeeeeteeueuuuu«auuuuuuuuc
C C
CACCCCCCCC DATA BLOCK BLOCK 0200
C C

DATA EPSEPSSJINJOUT/l.0D-30.0.(XX)1D+0.5.7/
C C
can“«euueuar-ar-uscuMun-«eeeeuueuuuuuueuuuuuuuc
C C
CICCCCCCCC INITIALIZATION BLOCK BLOCK 0400
C C

WRTTE(*.")’ENTER TIE NAME OF TIE DATA FILE’

READ(*.’(A40)’) DFILE

OPEN(8.FILE=DFILE)

WRTTE("'."')’ENTER TIE NAME OF TIE OUTPUT FILE’
C

READ(*.’(A40)’) OUTFIL

OPEN(7.FILE=OUTFIL)
C C
Car-unusua-aecu“ueuear-escenetetanus-«unusua-utens-Mute
C C
CPCCCCCCCC PROCESS BLOCK BLOCK 0500
C C
C -- START INPUT
C BLOCK 1

WRITE(7.*)’BEGIN LISTING INPUT QUANTITES’
200 READ(8,") N.NP.NT.ITMAX.MODL.IPRINT

WRITE(7,*)

WRITE(7.*)’BLOCK 1’

WRITE(7.")’N = NO. DATA POINTS. NP = NO. PARAMETERS’
WRITE(7.*)’NT = NO. OF INDEPENDENT VARIABLES’
WRITE(7.")’I'TMAX = MAXIMUM NO. OF ITERATIONS’
WRITE(7.*)’MODEL = MODEL NUMBER. IF SEVERAL MODELS IN SUBROUTINES:
] MODEL AND SENS’

WRITE(7,")’IPRINT = 1 FOR USUAL PRINTOUTS. 0 FOR LESS’
WRlTE(7.*)

IF(N.LE.0) TIEN

STOP
END IF
*.’(/.9X.”N".8X,”NP”.8X,"NT”.5X.”ITMAX”,5X,

+”MODEL”,4X.”IPRINT”)’)

WRITE(*,’(7IlO)’) N.NP.NT.ITMAX.MODL.IPRINT
WRITE(7.’(/.9X.”N".8X.”NP”.8X."NT”.5X.”ITMAX”.5X.
+”MODEL”.4X.”IPRINT”)’)

WRITE(7.’(7I]O)’) N.NP.NT.TTMAX.MODL.IPRINT

IOPT=0

C -- IF IOPT=O TIEN ON TIE 2ND AND SUCCEEDING ST ACKED CASES. TIE DATA

93

IS
C .. NOT REPRINTED.
C -- IF IPRINT=1. EXTRA PRINT OUT OF ETA. RESIDUALS B(]).... ARE GIVEN.
C BLOCK 2

WRITE(7."')

WRITE(7."')’BLOCK 2’

WRITE(7."')’B(1).B(2)....B(NP) ARE INITIAL PARAMETER ESTIMATES’

WRITE(7.")

READ(8.*)(B(I).I=1.NP)

WW7.'(10X."B(".I1.”) = ”.FlG.5)') (13(DJ=1.NP)

DO 150 J 1=2.5
BSUI) = O
150 CONTINUE
C
IF(IOPT 1.5.0) TIEN
C BLOCK 3
WRITE(7.")
WRITE(7.*)’BLOCK 3'
WRTTE(7."‘)’J = DATA POINT INDEX, Y0) = MEASURED VALUE'
WRITE(7.*)'SIGMA(I) = STANDARD DEVIATION OF YO)’
WRITE(7.*)'T(J.I) =- FIRST INDEPENDENT VARIABLE’
WRITECN’)
WRITE(7.'(/.9X."J".6X."Y(I)",3X."SIGMAU)".6X."T(J.l)"
+.6X."T(J.2)")')
D0 10 IZ=1.N
READ(8.")J.Y(J).SIGZ(I). (J.K1').KT=1.N'I)
W7o'0103F105Y) J.Y(J).SIG2(J).(T(J.KT).KT=1.NT)
8162(1) = SIGZ(J)*SIGZ(I)
TO CONTINUE
END IF
C
313 D0 2 IP=1.NP
DO 2 KP=1.NP
”(KPJP) ‘-" 0
P(KP.IP) = 0
CONTINUE
W7.'(/.SX."P(1.KP)".9X."P(2.KP)".9X."P(3.KP)”.9X.
+"P(4.KP)".9X."P(5.KP)")')
DO 6 IP=LNP

READ(8.")(PS(1P.KP).KP=1.NP)

WRITE(7.’(5D16.5)’) (PS(IP.KP).KP=],NP)
CONTINUE
BLOCK 4
DO 88 IP=1.NP
88 PS(IP,IP)=B(IP)“B(IP)*].OD+6

READ(8.*)EXTRA

C EXTRA=0 FOR NO EXTRA INPUT WHICH COULD BE FOR CONST ANTS

0°‘000000N

94

C IN MODELS
C =1 FOR ONE INPUT, NAMELY: EXTRA(1), ETC.
WRITE(7,")
WRITE(7,*)’BLOCK 4’
WRI'I'E(7.“)’IEXTRA = NO. OF EXTRAO) PARAMETERS. 0 IF NONE’
WRITE(7,"')
WRITE(7,’(10X,”IEXTRA = ”,IlO)’)IEXTRA
IF(IEXTRA .LT. 1) GOTO 21

WRIT'E(7."')
WRITE(7,*)’BLOCK 5’
WRI'I’E(7,*)’EXTRA(1),... ARE EXTRA CONSTANTS USED AS DESIRED’
WRITE(7,*)
READ(8,*)(EXTRA(IE),IE=1.IEXTRA)
W7.’("EX'1'RA(”.I2.") = ”.F16.S)’) (IBEXTRAW)E=1
1.IEXTRA)
21 CONTINUE
C
Cu-ADDBLANKCARDAFI'ERLASTINPUTCARD
C mEND INPUT
WRITE(7.*)’END INPUT QUANTITIES - - BEGIN OUTPUT CALCULATIONS’
WRITE(7.")
WRITE(7,*)’SY = SUM OF SQUARES FOR PRESENT PARAMETER VALUES’
WRITE(7.*)’SYP = SUM OF SQUARES FOR GAUSS PARAMETER VALUES.
SHOULD
I BE SMALLER THAN SY’
WRITE(7.*)’ SYP DECREASES TOWARD A POSITIVE CONSTANT’
WRITE(7,*)’G = MEASURE OF THE SLOPE, SHOULD BECOME SMALLER AS
llT’ERAT’IONS PROCEED’
WRITE(7.*)’ G SHOULD APPROACH ZERO AT CONVERGENCE’
WRITE(7,")’H = FRACTION OF THE GAUSS STEP. AS GIVEN BY THE
1BOX-KANEMASU METHOD’
WRITE(7,*)
WRITE(7.*)
DO 18 1L=1,NP
BSGLFBGL)
CG(IL) = O
18 CONTINUE
DO 19 IP==1.NP
XTY(IP)=0.0D+0
DO 19 KP=1.NP
HKPJP) = PS(KP.IP)
XTX(IP,KP)=0.0D+O
l9 CONTINUE
I = 0
MAX = O
C

99 MAX=MAX+1
C -- START BASIC LOOP GIVES B(I) AND SY
C

20

29

30

50

41

95

SY = 0.0D+0
DO 100 I3=1.N
1 = 13
CALL SENS
CALL MODEL
RISD = Y(l)-ETA
SY = SY + R18D‘RISD/SIG2(I)
SUM = 0.0D+0
DO 20 K=1.NP
XTY(K)=X'TY(I()+Z(K)*RISD/SIG2(I)
DO 20 L=1.NP
SUM = SUM + Z(L)"P(K.L)*Z(K)
XTX(K.L)= XTX(K.L) + Z(L)"Z(K)/SIGZ(D
CONTINUE
DELTA = 81020) + SUM
DO 29 JJ=1.NP
A0!) = 0.0D+O
CONTINUE
DO 30 JA=1.NP
DO 30 KA=1,NP
A(JA) = A(JA) + Z(KA)*P(JA,KA)
CONTINUE
C8 = 0.0D+0
DO 40 IC=1.NP
CS = CS + Z(JC)*(B(JC)—BS(IC))
(ISUC) = CG(IC) + Z(JC)*RISD/SIGZ(I)
CONTINUE
C = Y(I) - CS - ETA
DO 50 IB=1.NP
B(IB) = B(IB) + (ACIB)*C)/DELTA
CONTINUE
DO 41 ISV=1.NP
DO 41 JSV=1.NP
PSV(JSV.ISV) = P(JSV,ISV)
CONTINUE
DO 52 IV=1.NP
DO 52 IU=IV.NP
SUMP = 0.0D+O
DO 51 KP=1,NP
DO 51 JP=1.NP
IF(KP-IV.EQ.0.0R.JP-IU.EQ.0) GOTO 51
PSQ] = PSV(KPJP)*PSV(TU,IV)
PSQZ = PSV(IU,KP)"PSV(IV,JP)
PSQ = PSQl - PSQ2
IF(DABS(PSQ1)+DABS(PSQ2).LT.1.D-15) THEN
RP = PSQ "' 1.D15
ELSE
RP = PSQ / (DABS(PSQ1)+DABS(PSQ2))
END IF

RP = ABS(RP)
RPP = RP - 1.0D-12
IF(RPPLE.0.0D+O) THEN
PSQ = 0.0D+0
END IF
SUMP = SUMP + Z(JP)"Z(KP)“PSQ
51 CONTINUE
P(IU.IV) = (PSV(IU.IV)“SIG2(I)+SUMP)/DELTA
52 CONTINUE
DO 53 IV=2.NP
IVM = IV - 1
DO 53 IU = ”W
P(1UJV)= POVJU)
53 CONTINUE
IF(IPRINT.GT.O) THEN
IF(IEQ.1) THEN
WRITE(7.*)
WRITE(7,*)’SEQUENTIAL ESTIMATES OF THE PARAMETERS GIVEN BELOW’
WRITE(7,’(/l,3X,”I",6X,”ETA".5X.”RESIDUALS”,7X.
1"B(1)".8X."B(2)".8X."B(3)”.8X.”B(4)”)’)
END IF
WRITE(7,’(I4.6E12.5)’)IETA.RISD.(B(JC)JC=1.NP)
END IF
100 CONTINUE
C -- END BASIC LOOP. GIVES B(I) AND SY
C -- START BOX-KANEMASU MODIFICATION
C
C START BOX-KANEMASU MODIFICATION
IF(MAX-1)104,104.103
103 SS=SY/2.0D+0
IF(SS-SYP)104,104.105
105 D0 210 IBS=1.NP
B(IBS)= BSV(IBS)
210 CONTINUE
WRTTE(IOUT.212)
212 FORMAT(7X.'USE BSV(IBS)’)
GOT O 211
104 CONTINUE
DO 102 IBklNP
BSS(IBS)= BS(IBS)
102 CONTINUE
ALPHA= 2.0D+O
AA= 1.1D+0
110 ALPHA:- ALPHA/2.0M
DO 116 IBS=1.NP
BS(IBS)= BSS(IBS) + ALPHA“( B(IBS)-BSS(IBS) )
BSV(IBS)= BS(IBS)
116 CONTINUE
INDEX=0

V

G: 0.0D+0
DO 115 R134?
DELB= BS(IP)—BSS(TP)
G= G + DELB’CGCTP)
RATIO: DELB/( BSSCIP)+EPS )
RATIO= ABS(RATIO)
IF(RATIO-EPSS)1 13,113,114
113 INDEX: INDEX-H
WRITE(IOUT.314)
314 FORMAT(7X.’MAX’.8X.’NP’.5X,’INDEX’.8X.’IP’)
WRITE(7.'(7110)') MAX.NP.INDEX.IP
114 CONTINUE
C WRITE(7,122) I.Y(I).ETA.RISD.Z(IP),XYP.DELB,SIG2(I)
115 CONTINUE
SYP= 0.0D-I-0
DO 117 I3=1,N
I=13
CALL MODEL
RISD= Y(I)-ETA
SYP= SYP + RISD‘RISD/SIGZG)
117 CONTINUE
IF(NP-INDEX)106.106,107
106 H=1.0D+0
GOT O 132
107 CONTINUE
SYN: SYP*0.999D+0
IF(SYN-SY)112.112.111
111 IF(ALPHA-0.01D+0)109,109,110
109 WRITE(7.108) ALPHA,SYP,SY
108 FORMAT(3X.’ALPHA TOO SMALLAIPHA=’.F12.6.2X,’SYk’ElSfilX.
1'SY’,E15.6)
WRITE(7,1001)
1(X)1 FORMAT(BX.’R1)’,10X,’Z(2)’,10X,’Z(3)’.10X,’Z(4)’,10X.’Z(5)’)
1W FORMAT(6E13.4)
DO 1W3 1:1,N
CALL SENS
WRITE(7,1(X)2) (Z(IBB).IBB=1.NP)
1W3 CONTINUE
GOT O 1000
112 CONTINUE
SKSUM= SY - ALPHA*G*( 2.0D+0-l.0D+0/AA )
IF(SYP-SKSUM)131,131,130
130 H: ALPHA “ ALPHA‘GK SYP-SY+2.0D+0*ALPHA“G )
GOTO 132
131 CONTINUE
H= ALPHA*AA
132 CONTINUE
DO 118 IBN= LN?
B(IBN)g BSSGBN) + H " ( B(IBNIBSSGBN) )

 

 

98

118 CONTINUE
211 CONTINUE
WRITEOOUTJZI)
WRITE(*,I21)
121 FORMAT(SX,'MAX’,10X,’H’,13X,’G'.12X.
1’SY’.11X,’SYP’)
WRIT'E(7.122) MAXJ'I,G,SY,SYP
WRITE(*,122) MAXJ‘I.G.SY.SYP
122 FORMAT(IS,1F13.6,4E14.6)
WITIOXHBC'JL") = "$16.6” (1.3(DJ=1.NP)
W'.'(10x."3(".11.") '3 H.1316.6)') 03(0J319NP)
C END BOX-KANEMASU MODIFICATION
WRTTE(7,’(/.5X.”P(1,KP)”,9X,”P(2,KP)".9X.”P(3,KP)",9X,
1' 'P(4.KP)’ ’.9X.’ 'P(5.KP)' 0')
DO 206 IP=1.NP
WRIT'E(7.207) (POPJCP),KP=1,NP)
206 CONTINUE
207 FORMAT(5D15.7)
WRTTE(7,135)
135 FORMAT(SX,’CORRELATION MATRIX’)
DO 136 IR=1.NP
DO 136 IR2:=1.IR
AR: P(IR.1R) " P(IR2.IR2)
RURJRZF HRJRZVSQRTIAR)
136 CONTINUE
DO 137 IR=1.NP
WRITE(7.’(5E15.7)’) (R(IR,III),III=1,IR)
137 CONTINUE
DO 126 IPS=1.NP
PS(IPS.IPS)= (1-054'7) 1' MIPS)
126 CONTINUE
WRTI'E(7.*)’XTX(1.K).K=1.NP’
DO 220 K=1.NP
220 WRTTE(7.’(5E15.7)’XXTX(K.III),III=1,NP)
' WRITE(7,*)'XTY(I).I=1,NP'
WRTTE(7,’(5E15.7)’XXTY (1),I=1 ,NP)
127 FORMAT(3X,’IPS='.I4.3X.'PS(IPS.IPS)=',D15.8)
DO 119 IP=1.NP
XTY(IP)=0.0D+0
DO 119 KP=1.NP
1’01"ka PSGPKP)
XTX(IP.KP)=0.0D+0
119 CONTINUE
DO 120 IP=1.NP
BSGPF B(IP)
CG(IP)= 0.0D+0
120 CONTINUE
WRTIE(7,314)
WRITE(7,’(7110,4F10.4)’) MAXNPJNDEXJP

IF(NP-INDEX)101,101,123

123 CONTINUE
M=ITMAX
IF(MAX-M)99,99,IOI

101 CONTINUE
IF(IPRINT)I33,133.134

133 IPRINT =IPRINT +1
GOT O 99

134 CONTINUE

C

1(XX) CONTINUE
CLOSE(IIN)
CLOSEOOUT)

C

Ctttttttttttttttttttttt##ttttitttttttttttfittitttttttttttttttttttttc

CECCCCCCCC ERROR MESSAGES
E
(junta-nutnuutuunuuuuuuunuuuuuuuututu-u«C
CFCCCCCCCC FORMAT STATEMENTS
S

CtttttttIIttitttttttttttfittttitttlfitttitttit*tttttitttttittttitttttc

C

O

BLOCKOWO

GOO

BLOCK 9000

000

STOP
END
SUBROUTINE MODEL
C THIS SUBROUTINE IS FOR CALCULATING ETA, THE MODEL VALUE
IMPLICIT REAL‘8 (A-H,O-Z)
DIMENSION T(3500.5).Y(35(X)),SIG2(35(X1),B(5)Z(5),
+A(5).BS(5),VINV(5.5).EXTRA(20)
DIMENSION P(5.5).PS(5,5)

C INT'I-IISPROGRAM,THETA=ETA ANDETA=BETA
DOUBLE PRECISION KSL, Q. L. PI, ALPHSL
DOUBLE PRECISION ERFC, ZBRENT
DOUBLE PRECISION ALPHSR. PISR, X. XX. X2. XZALPH.
+ XALPI-I, XXALPH, BETA. BETAS. BETA2. T'HETA. ETA

COMMON $162,122.88,I,ETA,PS.P.B.A.Y.MODL.VINV.NP
+EXT'RA
COMMON/MOD/AA.'I‘L
COWON/PROP/KSL. Q, ALPHSL, L, PI
C
EXTERNAL ERFC. ZBRENT
C
KSL = 1.000

0000

100

ALPHSL = 1.0D0
L = BS(I)

Q = -1.0D0

BETA = T(I,1)

P1 = DACOS(-1.D0)

ALPHSR = ALPHSL“0.5DO
PISR = (PI“O.SDO)/2.D0

X 18 THE CALCULATED VALUE FOR LAMBDA
ZBRENrisamotfimmgsubmufineusedmsolvememscmdentalequafionforthe
freezing from location. See Numetical Recipes by Press et al., Cambridge University
Ptess. New York. New Yank. 1986.

X = ZBRENT(1.0D-4.2.D0,0.001D0)
XX = X*X
X2 = 2.DO*X
X2ALPH = X2‘ALPI-ISR
XALPH = X‘ALPHSR
XXALPH = XX*ALPIISL
CALCULATION OF DIMENSIONLESS TEMPERATURES
BETAS = BETA‘BETA
BETA2 = 2.D0*BET‘A
CALCULATION OF FROZEN PORTION TEMPERATURE
IF(BETA .LT. X) THEN
THE’I‘A = l-Q*(EXP(-BETAS)/BETA2 - EXP(-XX)/X2
-PISR‘(ERFC(BETA) - ERFC(X)»
ELSE
CALCULATION OF UNFROZEN PORTION TEMPERATURE
IF(BETA .GT. X) TI-IEN
THETA = (EXP(-ALPHSL*BETAS)/(BETA2*ALPHSR)
- PISR‘ERFC(ALPHSR*BETA))/(EXP(-XXALPII)/X2ALPH
-PISR*ERFC(XALPH))
ELSE
TEMPERATURE AT THE INTERFACE, DETERMINED FROM B.C.
TI-IETA = l.D0
ENDIF
ENDIF
ETA = THETA
RETURN
END
CALCULATION OF LAMBDA FROM FUNCTION ZBRENT
DOUBLE PRECISION FUNCTION ZBRENT(X1, X2, TOL)

PARAMETER(1TMAX=1(X), EPS= 3.0E-8)
DOUBLE PRECISION A. B, C. D. E, FA. FB. FC
DOUBLE PRECISION TOLI, TOL. X1, X2, XM
DOUBLE PRECISION P, Q. R. S, FUNCL

EXTERNAL FUNCL

+

101

=X1
B=X2
FA=FUNCL(A)
FB=FUNCL(B)

IF(FB*FA .GT. 0.0D0) PAUSE ’ROOT MUST BE BRACKETED FOR
ZBRENT.’
FC=FB
DO 15 ITER=1,TTMAX
IF(FB‘FC .GT. 0.0D0) THEN
GA
FCkFA
D=B-A
E=D
ENDIF
1F(ABS(FC) .LT. ABS(FB)) THEN
A=B
B=C
C=A
FA=FB
FB=FC
FC=FA
ENDIF

TOL1=2.0DO*EPS*ABS(B)+0.5DO*TOL
XM=0.5D0"(C—B)
IF(ABS(XM) .LE. TOLI .OR. FB .EQ. 0.0D0)TIEN
ZBRENT=B
RETURN
ENDIF

1F(ABS(E) .GE. TOLI .AND. ABS(FA) .GT. ABS(FB)) THEN
S=FBIFA
IF(A .EQ. C)THEN
P=2.0D0"XM"'S
Q=1.0D0 - S
ELSE
Q=FAIFC
R=FB/FC
P=S*(2.0D0"XM*Q*(Q-R) - (B-A)*(R-1.0D0))
Q=(Q-1.0IX))"'(R-I.0D0)*(S-1.0D0)
ENDIF
IF(P .GT. 0) Q = —Q
P=ABS(P)
IF(2.0D0*P .LT. MIN(3.0DO*XM*Q-ABS(T‘OL1"Q),ABS(E*Q)))TIIEN
ED
D=P/Q
ELSE
D=XM

OD

102

E=D
ENDIF
ELSE
D=XM
E=D
ENDIF
A=B
FA=FB
1F(ABS(D) .GT. TOL1)T'HEN
B=B+D
ELSE
=B+SIGNCTOL1,XM)
ENDIF
FB=FUNCL(B)
(DNTINUE

PAUSE ’ZBRENT EXCEEDING MAXIMUM HERAT'IONS.’
ZBRENT=B

RETURN

END

DOUBLE PRECISION FUNCTION ERFC(X)

DOUBLE PRECISION A1, A2, A3, A4, A5, P. T, X
A1=0.254829592D0
A2=—0.284496736D0
A3=1.421413741D0
A4=-1.453152027D0
A5=1.061405429D0
P=0.3275911D0
T=1.0D0/(1 .0D0+P*X)
ERFC=(AI *T+A2"T“2.0D0+A3"T"*3.0D0+A4*T*"'4.0D0+A5 *T“5.0D0)
*EXP(-X**2.0D0)
RETURN
END

DOUBLE PRECISION FUNCTION FUNCL(X)

DOUBLE PRECISION KSL. Q. ALPHSL, L, X, PI, ERFC
DOUBLE PRECISION EXPX. EXPXASL. XX2. RATIO
COMMON/PROP/KSL. Q. ALPHSL. L. PI

EXTERNAL ERFC

EXPX=EXP(-X"X)
EXPXASIzEXH-X‘X‘ALPHSL)
XX2=X*X*2.0D0
RATIO=EXPXASIJ(XX2*ALPHSL“0.5D0)

C

103

FUNCL-"KSL‘Q‘EXPX/XX2 - RATIO/(RATIO—(PI”0.5D0/2)
l"ERFC(ALPHSL'“"'0.5D0"'X)) -L"X

RETURN
END

SUBROUTINE SENS

C THIS SUBROUTINE IS FOR CALCULATING TI'IE SENSITIVITY COEFFICIENTS

C

0

00000

IMPLICIT REAL‘8 (A-H,O-Z)

DIMENSION T(3500.5),Y(35(1)),SIG2(35(X)).B(5).
+Z(5).A(5).BS(5).V1NV(5.5).EXTRA(20)
DIMENSION P(5.5).PS(5.5)

DOUBLE PRECISION KSL, Q. ALPHSL, L. PI
DOUBLE PRECISION ERFC, ZBRENT
DOUBLE PRECISION ALPHSR. PISR, X. XX, X2, X2ALPH.
XALPH, XXALPH, BETA, BETAS, BET A2, BOTTOM, TOP.
TOP2, DER], DER2A, DER2, DER3. DER4. DERS.
DER6. DER7, DER8, SENSEI

COMMON SIGZ.T.Z.BSJETAPSPEAYMODLNINVNP
+.EXTRA
COMMON/MOD/AA.TL

COMMON/PROP/KSL. Q. ALPHSL. L. P1

EXTERNAL ERFC, ZBRENT
KSL = 1.0D0
ALPHSL = 1.0D0
L = BS(I)
Q = -I.0D0
P1 = DACOS(-1.D0)
BETA = T(I,1)

VARIABLES DECLARED
ALPHSR = ALPHSL”O.5D0
PISR = (PI“0.5D0)/2.D0.

X IS THE CALCULATED VALUE FOR LAMBDA
ZBRENTisamotfindingsubrouflmnsedmsolvetheuanscendemalequafionforthe
freezing from location. See Numerical Recipes by Press et al., Cambridge University
Press, New York. New York. 1986.

X = ZBRENT(1.0D-4.2.D0,0.001D0)
XX = X‘X
X2 = 2.DO*X
X2ALPH = X2'ALPIISR
XALPI-I = X‘ALPHSR
XXALPH = XX*ALPI-ISL
BETAS = BETA*BETA
BETA2 = 2.0DO*BETA
BOTTOM = (EXP(-XXALPII))/X2ALPH - PISR*ERFC(XALPH)
TOP = (EXP(-XXALPH))/(2.0D0*XX"ALPHSR)
TOP2 = (EXP(-ALPHSL*BETAS)I(BETA2*ALPIISR)) - PISR

000000000000

00000

+

+++++

+++++

++

++++++

104
‘ERFC(BETA’ALPHSR)

SENSEI = TIE SENSITIVITY COEFFICENT FOR TIE TEMPERATURE
WITH RESPECT TO L

SENSE2 = TIE SENSITIVITY COEFFICENT FOR TIE TEMPERATURE
WITH RESPECT TO Q. NOT USED

SENSE3 = TIE SENSITIVITY COEFFICENT FOR TIE TEMPERATURE
WITH RESPECT TO KSL

SENSE/I = TIE SENSITIVITY COEFFICENT FOR TIE TEMPERATURE
WITH RESPECT TO ALPHASL

CALCULATION OF TIE SENSITIVITY COEFFICENTS

DERI = Q*(-EXP(-XX)/(2.0D0“XX))

DERZA = KSL*Q*((-XX‘EXP(-XX) - EXP(-XX))/(X*X*X)) -
(((-XX*EXP(-XXALPH) - EXP(-XXALPH))/(X*X*X*
ALPHSR))*BOTT OM - TOP‘(((-2.0D0*XX*ALPHSL**
(3.0D0/2.0D0)*EXP(-XXALPH) - ALPHSR*EXP
(-XXALPH))/(2.0DO"XX*ALPHSL)) + ALPHSR‘EXP
(-XXALPH))/(BOTTOM*BOT'TOM)) - L

DER2 = -DER2A“""(-I.ODO)

DER3 = -

DER4 = KSL“(EXP(-XX)I(2.0D0")O())

DEF-5 = Q*(EXP(-XX)/(2.0D0*XX))

DER6 = -(((-2.0D0*XX*ALPHSR*EXP(-XXALPH) - ALPHSL

""(—0.5D0)*EXP(-XXALPH))/(4.0D0"XX*ALPHSL))

‘BOT'TOM - TOP“((-2.0D0*XX*ALPIISR*
EXP(-XXALPH) - ALPHSL**(-0.5D0)*EXP
(-XXALPH))/(4.0D0"X*ALPHSL) + 0.5D0*ALPHSL
“(-0.5D0)*X*EXP(-XXALPH)))/(B OTTOM*BOTTOM)

DER7 = -TOP2*((-2.0D0“XX*(ALPHSL"'*(3.0D0/2.0D0))*EXP
(-XXALPH)) - (ALPHSR‘EXP(-XXALPH))/(2.0D0*XX*ALPHSL)
+ ALPHSR*EXP(-XXALPH))/(BOTTOM*BOTT OM)

DER8 =-(((-2.0D0*BETAS*ALPHSR‘EXP(-ALPHSL*BETAS) -
((ALPHSL“*(—O.5D0))"EXP(-ALPHSL‘BETAS)))/(4.0D0
*BETA*ALPHSL) + (0.5D0‘(ALPHSL“(—0.5D0))"‘BETA*
EXP(-ALPHSL"BET AS)))"BOTT OM - TOP2“‘(((-2.0D0"XX"I
ALPHSR*EXP(-XXALPH)) - ((ALPHSL“(-0.5D0))*EXP
(-XXALPH)))/(4.0D0*X*ALPHSL) + (0.5D0*"ALPHSL"
(-0.5D0))*X*EXP(-XXALPH)))/(BOTT OM‘BOTT OM)

CALCULATION OF FROZEN PORTION SENSITIVITY COEFFICENTS
IF(BET A .LT. X) TIEN
SENSEI 3' DERI‘DERZ‘DER3
SENSE2 = -(EXP(-BETAS)IBET A2 - EXP(-XX)/X2 - PISR

+ *(ERPQBETA) - ERFC(X)» - Q*(DER1*
+ DER2*DER4)

SENSEI = DERI‘DERZ‘DERS
SENSE2 = DERI'DERZ‘DER6
ELSE

000 0 000 00

00

105

CALCULATION OF UNFROZEN PORTION SENSITIVITY
COEFFICENTS
IF(BETA .GT. X) TIEN
SENSE] = DER7‘DER2‘DER3
SENSE2 = DER7‘DER2‘DER4
SENSE] = DER7‘DER2‘DER5
SENSE2 = DER8 + (DER7‘DER2‘DER6)
ELSE
SENSITIVITY COEFFICIENTS AT TIE INTERFACE
SENSE] = 0.0D0
SENSE2 a 0.0D0
SEN SE3 = 0.0D0
SENSE4 = 0.000
ENDIF
ENDIF
2(1) = SENSE]
2(2) = SENSE2

RETURN
END

SUBROUTINES MODEL AND SENSE MODIFIED TO
INCLUDE PRIOR INFORMATION

To include prior information in the subroutine MODEL: when the index equals 1, the

calculated value of theta is set equal to the estimated value of the parameter. The measured value

is set equal to the actual value of the parameter. obtained from prior information

To include prior information in the subroutine SENSE: when the index equals 1. the

sensitivity coefficient is set equal to 1.0.

SUBROUTINE MODEL

C THIS SUBROUTINE IS FOR CALCULATING ETA. TIE MODEL VALUE

IMPLICIT REAL*8 (A-II.O-Z)
DIMENSION T(35m.5).Y(35W).SIG2(35(XI).B(5).Z(5).

+A(5).BS(5).VINV(5.5).EXTRA(20)

DINENSION P(5.5)J’S(5.5)

C WRITTEN BY JAMES V. BECK

C

REVISED BY LESLE A. SCOTT
DOUBLE PRECISION KSL. Q. L. PI. ALPHSL

0

00000

106

DOUBLE PRECISION ERFC. ZBRENT
DOUBLE PRECISION ALPHSR. PISR. X. XX. X2. X2ALPH.
+ XALPH. XXALPH. BETA. BET AS. BET A2. TIETA. ETA

COMMON SIGZ.T.Z.BSJEIAPSPB.A.Y.MODL.VINV.NP
+.EXTRA
COMMON/MOD/AAJ'L
COMMON/PROP/KSL. Q. ALPHSL. L. Pl

EXTERNAL ERFC. ZBRENT

KSL = 1.0D0
ALPHSL = 1.0D0
L = BS(1)

Q = -].0D0

BETA = T(T.])

P1 = DACOS(-1.DO)

ALPHSR = ALPHSL“0.5D0
PISR = (PI“0.5D0)/2.D0

X IS THE CALCULATED VALUE FOR LAMBDA
ZBRENT is a root finding subroutine used to solve the transcendental equation for the
freezing front location. See Numerical Recipes by Press et al.. Cambridge University
Press. New York. New York. 1986.
X = ZBRENT(1.0D-4,2.D0.0.00IDO)
XX = X*X
X2 = 2.DO*X
X2ALPH = X2*ALPHSR
XALPH = X*ALPHSR
XXALPH = XX *ALPHSL
CALCULATION OF DIMENSIONLESS TEMPERATURES
BETAS a BETA'BETA
BETA2 = 2.D0*BET‘A
IF(I .EQ. 1)THEN
T'HETA = 88(1)
ELSE
CALCULATION OF FROZEN PORTION TEMPERATURE
IF(BETA .LT. X) THEN
THETA = I-Q*(EXP(-BETAS)/BETA2 - EXP(-XX)/X2
+ -PISR'(FRFC(BETA) - ERFC(X)»
ELSE
CALCULATION OF UNFROEN PORTION TEMPERATURE
IF(BET‘A .GT. X) THEN
T'HETA = (EXP(-ALPHSL'BETAS)/(BETA2"ALPHSR)
+ - PISR'ERFC(ALPHSR‘BETA))I(EXP(-XXALPH)/X2ALPH
+ -PISR*ERFC(XALPH))
ELSE
TEMPERATURE AT THE INTERFACE DETERMINED FROM B.C.

107

TIETA = I.D0
ENDIF
ENDIF
ENDIF
ETA = TIETA
RETURN
END
CALCULATION OF LAMBDA FROM FUNCTION ZBRENT
DOUBLE PRECISION FUNCTION ZBRENT(X1. X2. TOL)

PARAMETEROTMAX=1(X). EPS= 3.0E-8)
DOUBLE PRECISION A. B. C. D. E. FA. FB. FC
DOUBLE PRECISION TOL]. TOL. X1. X2. XM
DOUBLE PRECISION P. Q. R. S. FUNCL

EXTERNAL FUNCL
A=X1

B=X2
FA=FUNCL(A)
FB=FUNCL(B)

IF(FB’FA .GT. 0.0D0) PAUSE ’ROOT MUST BE BRACKETED FOR
ZBRENT.’
FC=FB
DO 15 ITER=1.ITMAX
IF(FB‘FC .GT. 0.0D0) TIEN
@A
FC=FA
D=B-A
E=D
ENDIF
1F(ABS(FC) .LT. ABS(FB)) TIEN
A=B
B=C
C=A
FA=FB
FB=FC
FC=FA
ENDIF

TOL1=2.0D0*EPS*ABS(B)+0.SDO"TOL
XM=0.5D0*(C-B)
1F(ABS(XM) .LE. TOL] .OR. FB .EQ. 0.0D0)TIEN
ZBRENT=B
RETURN
ENDIF

1F(ABS(E) .GE. TOL] .AND. ABS(FA) .GT. ABS(FB)) TIEN
S=FBIFA

.4.

108

IF(A .EQ. C)TIEN
P=2.0D0"XM"S
Q=].OD0 - S

ELSE
Q=FA/FC
R=FBIFC
P=S*(2.0D0*XM*Q*(Q-R) - (B-A)"(R-l.0D0))
Q=(Q-1.0D0)*(R-1.0D0)*(S-1.0D0)

ENDIF

IF(P .GT. 0) Q = -Q

kABS(P)

IF(2.0D0"P .LT. MTN(3.0DO“XM"Q-ABS(TOL1*Q).ABS(E*Q)))TIEN
E=D
D=P/Q

ELSE
D=XM
E=D

ENDIF
ELSE
D=XM
E=D
ENDIF
A=B
FA=FB
1F(ABS(D) .GT. TOL1)TIEN
B=B+D
ELSE
B=B+SIGN(TOLI .XM)
ENDIF
FB=FUNCL(B)
CONTINUE

PAUSE ’ZBRENT EXCEEDING MAXIMUM ITERATIONS.’
ZBRENT=B

RETURN

END

DOUBLE PRECISION FUNCTION ERFC(X)

DOUBLE PRECISION A]. A2. A3. A4. A5. P. T. X

A]=0.25482959ZDO

A2=-0.284496736D0

A3=].42]4]3741D0

Ah-l.453152027D0

A5=1.06l405429D0

P=0.327591 IDO

T=1.0D0/(1.0D0+P*X)
ERFC=(A]‘T+A2"T**2.0D0+A3‘T“3.0D0+A4*T“4.0D0+A5‘T"*5.0D0)

*EXP(-X“2.0D0)

00

109

RETURN
END

DOUBLE PRECISION FUNCTION FUNCL(X)

DOUBLE PRECISION KSL. Q. ALPHSL. L. X. PI. ERFC
DOUBLE PRECISION EXPX. EXPXASL. XX2. RATIO
COMMON/PROP/KSL. Q. ALPHSL. L. P]

EXTERNALERFC

EXPX=EXP(-X*X)
EXPXASIzEXP(-X*X*ALPHSL)
XX2=X*X*2.0D0
RATTO=EXPXASU(XX2*AIJ’HSL“0.SDO)
FUNCILa-KSL‘Q‘EXPX/XXZ - RATIO/(RATIO-(PI“0.5D0/2)
+ ‘ERFC(ALPHSL**0.5D0*X)) -L"'X
RETURN
END

SUBROUT’INE SENS

TIITS SUBROUTINE IS FOR CALCULATING TIE SENSITIVITY COEFFICENTS
IMPLICIT REAL'8 (A-H.O-Z)

DIMENSION T(35(X).5).Y(35M).SIG2(35(X)).B(5).
+Z(5).A(5).BS(5).VINV(5.5),EXTRA(20)

DIMENSION P(5.5).PS(5.5)

DOUBLE PRECISION KSL. Q. ALPHSL. L. P]
DOUBLE PRECISION ERFC. ZBRENT
DOUBLE PRECISION ALPHSR. PISR. X. XX. X2, X2ALPH.
+ XALPH. XXALPH. BETA. BETAS, BET A2. BOTTOM. TOP.
+ TOP2. DER]. DER2A, DER2. DER3. DER4. DERS.
+ DER6. DER7, DER8. SENSE]
COMMON S]G2.T.Z.BS.LE'I’AJ’SJ’B.A.Y.MODL,VII\IV.I‘II>
+.EX'I'RA
COMMON/MOD/AAJ‘L
COMMON/PROP/KSL. Q. ALPHSL. L. P]

EXTERNAL ERFC. ZBRENT
KSL = 1.0D0
ALPHSL = 1.0m
L = BS(I)
Q = -1.0D0
P1 = DACOS(-1.D0)
BETA = T(I.1)

VARIABLES DECLARED
ALPHSR = ALPHSL“‘0.5w
PISR = (PI“0.5D0)/2.D0

00000

0000000000

+++++

+++++

+

+

+

110

X 18 THE CALCULATED VALUE FOR LAMBDA
ZBRENTisamotfindingsubmufineusedmsolvemeuanscendemmequafionforme
freezing from location. See Numerical Recipes by Press et al., Cambridge University
Press. New York. New York. 1986.
X = ZBRENT(1.0D4,2.D0.0.00IDO)
XX = X‘X
X2 = 2.D0*X
X2ALPH = X2*ALPI-ISR
XALPH = X‘ALPHSR
XXALPH = XX‘ALPHSL
BETAS = BE‘I‘A*BET‘A
BETA2 = 2.0DO*BETA
BOTTOM = (EXP(-XXALPH))/X2ALPH - PISR‘ERFC(XALPH)
TOP = (EXH-XXALPH))/(2.0D0*XX*ALPHSR)
TOP2 = (EXP(-ALPHSL*BETAS)/(BETA2*ALPHSR)) - PISR
*ERFQBETA‘ALPHSR)

SENSE] = TIE SENSITIVITY COEFFICENT FOR TIE TEMPERATURE
WITH RESPECT TO L

SENSE2 = TIE SENSITIVITY COEFFICENT FOR TIE TEMPERATURE
WITH RESPECT TO Q. NOT USED

SENSE3 = TIE SENSITIVITY COEFFICENT FOR TIE TEMPERATURE
WITH RESPECT TO KSL

SENSE4 = TIE SENSITIVITY COEFFICENT FOR TIE TEMPERATURE
WITH RESPECT TO ALPHASL

CALCULATION OF TIE SENSITIVITY COEFFICENTS

DER] = Q‘(-EXP(-XX)/(2.0D0*XX))

DERZA = KSL‘Q*((-XX*BXP(-XX) - EXP(-XX))/(X*X"X)) -
(((-XX"EXP(-XXALPH) - EXP(-XXALPH))/(X*X*X*
ALPHSR))*BO’TTOM - TOP*(((-2.0D0*XX*ALPHSL**
(3.0D0/2.0D0)*EXP(-XXALPH) - ALPHSR*EXP
(~XXALPH))/(2.0D0"XX*ALPHSL)) + ALPHSR‘EXP
(~XXALPH))/(BOTTOM"'BOTTOM)) - L

DER2 = -DER2A“(-1.0D0)

DER3 :3 -X

DER4 = KSL'EXP(-XX)/(2.0DO"XX))

DERS = Q‘(EXP(-XX)I(2.0D0*XX))

DER6 = -(((-2.0D0*XX*ALPHSR"EXP(-XXALPH) - ALPHSL
I""(-0.5D0)"’EXP(-XXAIJ’H))/(4.0D0“'XX"ALPIISL))

*BOTTOM - TOP*((-2.0D0*XX*ALPHSR*

EXP(-XXALPH) - ALPHSL“(-0.5D0)"EXP
(-XXALPH))/(4.0D0*X*ALPHSL) + 0.5D0'ALPHSL
“(-0.5D0)"'X *EXP(-XXALPH)))/(BOTT OM‘BOTT’ OM)

DER7 = -TOP2"'((-2.0DO"XX"‘(ALPHSL“(3.0D0/2.0D0))"EXP
(-XXALPH)) - (ALPHSR‘EXP(-XXALPH))/(2.0DO*XX*ALPHSL)
+ ALPHSR’EXPGXXALPIDIKBOTT OM‘BOT'T OM)

DER8 =-(((-2.0D0"'BETAS*ALPHSR‘EXP(-ALPHSL*BE‘TAS) -
((ALPHSL“‘(-0.5D0))"‘EXP(-ALPHSL"BETAS)))/(4.0D0

0

000 0 000 00 00000

00

+++++

1]]

‘BE’TA‘ALPHSD + (0.5D0"'(ALPHSL""(-0.5D0))"'BE'I‘A"I
EXK-ALPHSL‘BETAS)))*BOTTOM - TOP2‘(((-2.0D0‘XX“
ALPHSR‘EXH-XXALPID) - ((ALPHSL“(-0.5D0))"EXP
(-XXALPH)))I(4.0D0*X*ALPHSL) + (0.5D0"'"'ALPIISL""‘l
(-O.5DO))"X*EXP(-XXALPII)))I(BOT'TOM*BOTT OM)

IF (I .EQ. 1) TIEN
SENSE] = 1.0D0
ELSE
CALCULATION OF FROZEN PORTION SENSITIVITY COEFFICENTS
IF(BETA .LT. X) TIEN
SENSE] = DER]“DER2*DER3
SENSE2 = -(EXP(-BETAS)/BETA2 - EXP(-XX)/X2 - PISR
+ ‘(ERFC(BETA) - ERFC(X)» - Q‘(DER1*
+ DER2‘DER4)
SENSE] = DER]*DER2"'DERS
SENSE2 = DER]*DER2"DER6
ELSE
CALCULATION OF UNI‘ROZEN PORTION SENSITIVITY
COEFFICENTS
IF(BETA .GT. X) TIEN
SENSE] = DER7‘DER2‘DER3
SENSE2 = DER7‘DER2‘DER4
SENSE] = DER7‘DER2‘DER5
SENSE2 = DER8 + (DER7‘DER2‘DER6)
ELSE
SENSITIVITY COEFFICIENTS AT TIE INTERFACE
SENSE] = 0.0D0
SENSE2 = 0.000
SENSE3 = 0.0D0
SENSE4 = 0.0D0
ENDIF
ENDIF
ENDIF
Z(1) = SENSE]
2(2) = SENSE2

RETURN
END

APPENDIX D

THE FORTRAN PROGRAM MODFOR

This program. MODPOR. is used to provide an input file for use with NLINA.FOR either

with or without random measurement enors.

PROGRAM MODEL

THIS PROGRAM IS DESIGNED TO DETERMINE TIE DIMENSIONLESS
TEMPERATURESAS FUNCTIONS OF POSITION AND TIME. OF BOTH
TIE FROZEN AND UNFROZEN REGIONS SURROUNDING A POINT
SOURCE IEAT SINK. WIIII RANDOM ERRORS

WRITTEN BY LESLE SCOTT

0000000

DOUBLE PRECISION KSL. Q, ALPHSL. L, PI
DOUBLE PRECISION ERFC, ZBRENT
DOUBLE PRECISION DETA, ALPHSR, PISR, X, XX, X2. X2ALPH.
+ XALPH. XXALPH. BETA. BET AS. BET A2. TIETA
DIMENSION DATA(2(XX)0)

COMMON/PROP/KSL, Q. ALPHSL. L. P]
COMMON/RAND/NT. S'TDDV

EXTERNAL ERFC. ZBRENT

OPEN(UNIT=10, FILE="TEMPS.DAT". ST ATUS="UNKNOWN")
KSL :8 ].D0

Q = -].ODO

ALPHSL = 1.000

L = ~1(X).DO

PI = DACOS(-].DO)

NT = 150

INCR = ]

DET A = 1.0D-2

ALPHSR = ALPHSL“0.5DO

]]2

113

PISR = (PI“0.SDO)/2.DO
C X 18 THE CALCULATED VALUE FOR LAMBDA
C ZBRENTisamotfindingsubmuflneusedmsolvetheuanscardemalequafionforme
C freezing from location. See Numerical Recipes by Press et al.. Cambridge University
C Press, New York. New York. 1986.
X = ZBRENT(1.0D-4.2.D0,0.001D0)
WRITE(10,*)"X=”,X
XX = X‘X
X2 = 2.D0*X
X2ALPH = XZ‘ALPHSR
XALPH = X‘ALPHSR
XXALPH = XX‘ALPHSL
C CALCULATION OF DIMENSIONLESS TEMPERATURES
C BETAISTI-IESAMEASETA
CALL RANDOM (DATA)
DO I = 2.N'I‘.INCR
BETA = (I-l)‘DETA
BETAS = BETA'BETA
BETA2 = 2.DO*BETA
C CALCULATION OF FROZEN PORTION TEMPERATURE
IF(BETA .LT. X) THEN
TI-IETA = 1-Q*(EXP(-BETAS)IBETA2 - EXP(-XX)/X2
+ -PISR‘(ERFC(BETA) - ERFC(X)»
ELSE
C CALCULATION OF UNFROZEN PORTION TEMPERATURE
IF(BETA .GT. X) THEN
THETA = (EXP(-ALPHSL*BETAS)/(BETA2"ALPHSR)

+ - PISR‘ERFC(ALPHSR*BETA))I(EXP(-XXALPII)/X2ALPII
+ -PISR‘ERFC(XALPII))
ELSE
C TEMPERATURE AT TIE INTERFACE, DETERMINED FROM B.C.
TIETA = ].D0
ENDIF
ENDIF

C ADDITION OF RANDOM ERRORS TO MEASUREMENT DATA
TIETA I: TIETA + DATAO-l)
WRIIE(]O,’(I lO.7F]O.5)')I-],TIETA.STDDV.BETA
ENDDO
STOP
END
C CALCULATION OF LAMBDA FROM FUNCTION ZBRENT
DOUBLE PRECISION FUNCTION ZBRENT(X], X2, TOL)

C
PARAMETERGTMAX=1m, EPS= 3.0E-8)
DOUBLE PRECISION A, B, C. D, E, FA, FB. FC
DOUBLE PRECISION TOL], TOL. X1, X2, XM
DOUBLE PRECISION P. Q, R, S, FUNCL

C

EXTERNAL FUNCL

+

114

A=X]
B=X2
FA=FUNCL(A)
FB=FUNCL(B)

IF(FB‘FA .GT. 0.0130) PAUSE 'ROOT MUST BE BRACKETED FOR
ZBRENT.‘ 4
FGFB
DO 15 ITER=].ITMAX
IF(FB‘FC .GT. 0.0130) TIEN
C=A
FC=FA
D=B -A
E=D
ENDIF
1F(ABS(FC) .LT. ABS(FB)) TIEN
A=B
B=C
C=A
FA=FB
FB=FC
FC=FA
ENDIF

TOL]=2.0D0‘EPS*ABS(B)+O.5D0*TOL
XM=O.5DO*(C-B)
1F(ABS(XM) .LE. TOL] .OR. FB .EQ. 0.0D0)TIEN
ZBRENT=B
RETURN
ENDIF

1F(ABS(E) .GE. TOLI .AND. ABS(FA) .GT. ABS(FB)) THEN
S=FBIFA
IF(A EQ. C)TIEN
P=2.0DO"XM*S
Q=l.0D0 - S
ELSE
Q=FAIFC
R=FBIFC
PP-S’(2.0DO*XM*Q"(Q-R) - (B-A)*(R-I.ODO))
Q=(Q-I.ODO)*(R-].0D0)"(S-].0D0)
ENDIF
IF(P .GT. 0) Q = -Q
P=ABS(P)
IF(2.0D0‘P .LT. MlN(3.0D0"XM*Q-ABS(TOL1*Q),ABS(E*Q)))TIEN
E=D
D=P/Q
ELSE
D=XM

115

E=D
ENDIF
ELSE
D=XM
E=D
ENDIF
A=B
FA=FB
IF(ABS(D) .GT. TOLDTIEN
B=B+D
ELSE
B=B+SIGN(TOL] .XM)
ENDIF
FB=FUNCL(B)
CONTINUE

PAUSE ’ZBRENT EXGEDING MAXIMUM ITERATIONS.’
ZBRENT=B

RETURN

END

DOUBLE PRECISION FUNCTION ERFC(X)

DOUBLE PRECISION A], A2, A3, A4, A5. P. T, X
A]=0.254829592m
A2=-O.284496736D0 '
A3=1.42141374]D0
A4=-].453]52027DO
A5=I.061405429D0
P=0.32759] IDO
T=l .0D0/( 1 .0D0+P"X)
ERFC=(A1*T+A2"T"*2.0D0+A3"T“3.0D0+A4"I‘"4.0D0+A5*T“5.0DO)
*EXP(-X“2.0D0)
RETURN
END

DOUBLE PRECISION FUNCTION FUNCL(X)

DOUBLE PRECISION KSL. Q. ALPHSL. L. X, PI, ERFC
DOUBLE PRECISION EXPX. EXPXASL. XX2. RATIO
COMMON/PROP/KSL. Q. ALPHSL. L. P]

EXTERNALERFC

EXPX=EXP(-X"X)

EXPXASL=EXP(-X*X*ALPHSL)

XX2=X*X*2.0DO
RATIO=EXPXASIJ(XX2*ALPHSL“O.5DO)
FUNCL=KSL"Q"EXPX/XX2 - RATIO/(RATIO-(PI“O.5D0/2)

0099 0

000 0

1]
12

116

+ *ERFC(ALPHSL“O.5D0‘X)) -L*X
RETURN
END

SUBROUTINE RANDOM (DATA)
COMMON/RAND/NT. STDDV
COMMON NDATNPI‘S
PARAMETER(PI=3.14159265.NPIS=4.NBm=1000.NDAT=NPTS+NBm)
PARAMETER(PI=3. 14159265 .NBIN=1000)
SEE Numerical Recipes by Press. Flanner'y, Teukolsky and Vetteriing.
Cambridge Press. 1986 about page 192
Modified by JV Beck. Michigan State University,
E-mail address: 22427jvb@itm.cl.msu.edu
DIMENSION DATA(2(D00)
CHARACI'ER*80 FOUT
NPI‘S = NT
IDUM IS SEED. SET TO ANY NEGATIVE NUMBER TO INITIALIZE OR
REINITIALIZE.
IDUM=-5
WRITE(*,*)’ Eruer the number of points ’
READ(*,")NPTS
WRI'IE(*,*)'ENTER THE SEED NUMBER (-)’
READ(*.*)IDUM
NDAT=NPTS+NBIN
WRITE(*,*)’ GIVE THE STANDARD DEVIATION’
READ(*,*)STDDV
WRITE(*,*)'Give the name of the output file'
READ(*.‘(A80)’)FOUT
OPEN(13. FILBFOUT)

RHON=0.0

RHOD=0.0

WRITE(*.*)' I RAND. NO.’

DO 500 IDUMI=1,1

DATA(1)=GASDEV(IDUM)"STDDV

WRITE(*,1(I))1.DATA(1)

WRITE(13.1(X))].DATA(1)

DO 1] I=2.NPTS
DATA(I)=GASDEV(IDUM)"STDDV
RHON=RHON+DATA(I-l)*DATA(I)
RHOD=RHOD+DATA(I)*DATA(I)

WRITE(*,]00)I,DATA(T)

WRITE(1 3,1CD)I.DATA(I)

CONTINUE
CONTINUE
RHO=RIION/RIIOD

CCC WRITE(",’(1X.A/)’) ’Descriptors of a gaussian distribution’

CALL MOMENT(DATA.I-LAVE.ADEV.SDEV.VAR.RHO)

SIDCONTINUE

C

117

WRII‘E(*.’(IX.T29.A.T42,A/)’) ’ Values of mantities’.’ ’
WRITE(*,*)' Values of quarrtities’
WRI'I‘E(*,’(1X,T29,A,T42.A/)’) ’ Sample ’,’Expected’
WRITE(*,’(IX,A,T25.2F12.4)’) ’Mean :’.AVE.0.0
WRITE(*,‘(1X.A.T25.2F12.4)’) ’Average Deviation :',ADEV,STDDV
WRITE(*.’(1X.A.T25.2F12.4)’) ’Standard Deviation :’.SDEV.STDDV
VARTH=STDDV*STDDV
WRIIE(*,’(1X.A.T25.2F12.4)’) ’Variance :’,VAR.VARTH
WRI'I'E(*.‘(1X.A,T25.F12.4)’)‘Est. Correlation Coef.'.RHO
WRITE(*,")’Average deviation comes from use of absolute values’

100 FORMAT(IIOJ" 10.6)

11

12

END
SUBROUTINE MOMENT(DATA,N.AVE.ADEV,SDEV,VAR,RHO)
DIMENSION DATA(2(X)00)
IF(NLE.1)PAUSE 'N must be at least 2’
S=O.
SD=O.
SN=0.
DO 11 J=1.N
S=S+DATA(J)
IF(J .EQ. 1)GO’TO ll
SN=SN+DATA(J)*DATA(J-l)
SD=SD+DATA(J)*DAT‘A(J)
CONTINUE
AVE=S/N
ADEV=0.
VAR=0.
DO 12 MN
S=DATA(J)-AVE
ADEV=ADEV+ABS(S)
P=S*S
VAR=VAR+P
CONTINUE
ADEV=ADEVIN
VAR=VAR/(N-l)
SDEV=SQRT(VAR)
RHO=SNISD
WRITE(*,*)‘SN SD RHO’SNSDRHO
RETURN
END
FUNCTION RAN 1(IDUM)
DIMENSION R(97)

RETURNS UNIFORMLY DISTRIBUTED NUMBERS BETWEEN 0 AND 1

PARAMETER (M1=2592CXIJA1=7141,IC1=54773.RM]=3.8580247E-6)
PARAMETER (M2=134456,IA2=8121.I(2=28411.RM2=7.4373773E-6)
PARAMETER (M3=243(XX).IA3=4561,IC3=51349)
DATA IFF DI
IF (IDUM.LT.0.0R.IFFEQ.O) TIEN

IFF=1

118

IX 1 =MOD(IC1 ~IDUM.M 1)
IX1=MOD(IA1"IX]+IC1.M])
IX2=MOD(IX] .M2)
IX]=MOD(IA1*IX1+IC1.M1)
IX3=MOD(IX1.M3)

DO 1 1 1:1,97
IX1=MOD(IA1"'IX1+IC1.M1)
IX2=MOD(IA2*IX2+IC2.M2)
R(J)=(FLOAT(IX1)+FLOAT(IX2)‘RM2)‘RM1

I 1 CONTINUE
IDUM=1
ENDE
IX]=MOD(IA]"IX1+IC1,M1)
IX2=MOD(IA2"IX2+IC2,M2)
IX3=MOD(IA3‘IX3+IC3.M3)
J=1+(97*IX3)/M3
E(J.GT.97.0R.J.LT. ])PAUSE
RAN ]=R(J)
R(I)=(FLOAT(IX 1)+FLOAT(IX2)*RM2)’RM1
C W*.‘)’J.R(J).RAN1’J.R(D.RAN1
RETURN

END
FUNCTION GASDEVCIDUM)
USES BOX-MULLER TRANSFORMATION FROM UNEORM DISTRIBUTION TO
NORMAL DISTRIBUTION WITI-I UNIT STANDARD DEVIATION
DATA [SET/OI
E (ISETEQD) TIEN
1 VI=2.*RAN1(IDUM)-1.
V2=2."RAN1(IDUM)-1.
R=V]“2+V2“2
E(R.GE.]..OR.R.EQ.O.)GO TO 1
FAC=SQRT(-2.’LOG(R)/R)
GSET=V]"'FAC
GASDEV=V2"'FAC
ISET =1
ELSE
GASDEV=GSET
ISE'T=0
ENDE
C WRIIE(",‘)’IDUM.GASDEV’.IDUM.GASDEV
RETURN
END

00

119
This file represents a sample output file from MOD.FOR. without prior information. to

be used as input for NLINA.FOR for estimation of dimensionless laterrtheat of fusion. L'. The
firstrowofmmrbersrepresemthemttnberofdatapoints,themrmberofparameterstobe
estimated, the munber of independent variables, the maximum munber of iterations to be
performed, the model number, and the usual printouts respectively. The second row represents
theinitialguessofl.’.whicbistobeestimated. Thefirstcolumnistheindex,thesecondcolumn
isthevaluesofthedimertsionlesstemperauues.drethirdisdrestandarddeviationofthe

meamrement errors, and the fourth column is the independent variable 11.

126.].1.](X).1,1
~150 0d0
1 48.02077 .orooo .0](XX)
2 23.04473 .0100) .02000
3 14.69203 .01000 .03000
4 10.52862 .OICXX) 04“”
5 8.03962 .01(XX) .05(XX)
6 6.38015 .01(XX) .06(XX)
7 5.21252 .orooo 07(0)
8 4.29658 .0101” .08(XX)
9 3.63488 .orooo .09(XX)
10 3.07533 .01000 .roooo
]1 2.60642 .01000 .11000
12 2.22666 010(1) .12000
13 1.93491 .01000 .13000
14 1.66083 01“!) .14000
15 1.42710 .01(XX) .15(XX)
16 1.21884 .010“) .16(XX)
]7 1.0550] .0](XX) .17GX)
18 .9666] 01% .]8(XX)
19 .88892 .01000 .190“) .
20 .82373 .OIIXX) .20000
21 .7591] 01“!) 21“!)
22 .69562 .OIIXX) 22(0)
23 .67533 .OIIXX) .2311!)
24 .63256 .01“ .2411)
25 .58378 .01(XX) .2511!)
26 .55002 .0101) .260!)
27 .52519 01“!) 27000
28 .49508 .OIGX) 28(0)
29 .48083 .0101) .291XX)
30 .4497] .01(XX) .30000
31 .41243 .0111!) 31m

32
33

35

37
38
39
40
4]
42
43

45
47
49

5]
52
53

55
56
57
58
59

61
62
63

65

67
68

70
7]

73
74
75
76

78

.39110
.38648
.35509

.3174]
.30356

.29075
.27858
.26758
.23265
.22042

.21148
.20164
.18726
.19548
.16875
.15592
.16928
.14798
.1487]
.13669
.14019
.14412
.12903
.12897
.11703
.09989
.11688
.10896
.09443
.08428
.09883

.08838
.00045
.08647
.08845

.07538

.07927
.04960
.05827
.04116
.07273
.05815
.03799

.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000
.01000

.32000
.33000

.35000
:37000
.38000
.39000

:4rooo

:57000
38000

.59000

.61000

120

81
82
83
84
85
86
87
88
89
90
91

93
94

98

]W
101
102
103
104
105
106
107
108
109
110
111
112
113
114
]15
116
117
118
119
120
121
122
123
124
125
126

02042

.05090
.05272
04.330

.02789
.03682

.04152
.05046
.03143
.0479]
.02869
.04304
.03614
.03 1 38

02306
.m666
.04897
.01989
.02609

.02754
.W568
.01499
.01834

.01425
.01647
.02279
.W334
.02032
.01718

.01142
.02530
0.1360

.02494

.01605

-.W265

.W633

- .W558

.(XXISO
.02777

«(X3227

.W57]
.02586
.W265
.W251
.W160

.01(XX)
.01(XX)
.01(XX)
.OIIXX)
.01(XX)
.01(XX)
.0111!)
.0111”
.01W0
.010W
.01(XX)
.01(XX)
.01(XX)
01“!)
.01000

.01W0
.010W
.01“
0““)
.010W
.01(XX)
.010W
.01W0
.01W0
.OICXX)
.010W
.01W0
.OICXX)
01(0)
.01”
.01(XX)
.01(XX)
.010W
.010W
.OIIXX)
.01(XX)
.01(XX)
.OICXX)
.01(XX)
.01(XX)
.011!”
01“!)
.01(XX)
.0111!)
.01(XX)
.0](XX)

.81(XX)
.82(XX)
.83CXX)

.85W0
.8611!)
.87CXX)
.8811”

.91W0
.9301)
.95(XX)

.97IXX)
.98(XX)

1 WOW
101(XX)
1.020W
1.03W0
1.040W
1.05W0
1.06(XX)
1.07W0
1.08W0
1.090W
1.10000
1.1]W0
1.12000
1.13“!)
1.]4W0
1.15111)
1.16W0
1.17W0
1.18“!)
1.19“!)
120000
121000
1.22“!)
1.23111)
124000
125000
1.26“!)

]21

122

This file represents a sample output file from MODFOR. with prior information. to be

used as input for NLINA.FOR for the estimation of L’.

]38.1.1.1W.1.]

-1W.(XX)
48.02888
23.03408
14.70647
10.54467
8.04848
6. 38618
5 .20235
4. 3141 3
3.62446
3.07378
2.625 17
2.25082
1.93 376
1.66656
1.43220
1.22878
1.04916
.94665
.8810]
.81855
.76560
.71388
.67053
.631 30
.593W
.55770
.52673
.49654
.47143
.44409
.42016
.401W
.3787]
.36145
.34176
. 32408
.30975
.29423
.27993
.26833
.25474
.2418]

10000

.W1W
.W1W
.W1W
.W1W

.W1W
.W1W
.W1W
.W1W

.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W
.W1W

.01(XX)
.OIGX)

.11(XX)
.12CXX)
.13CXX)
.14IXX)
.ISIXX)
.160W
.17(XX)

.18(XX)
.190W

.21W0
.2311”
25(0)
.27W0
.2811!)
29(0)
.3111!)
.32”
.3311!)
.35W0
.3711”
.381“)
.39(XX)

241000
.42000

This row contains the prior information of L’

44
45

47
49

51
52
53

55
56
57
58
59

61
62
63

67
68

70
71

73
74
75
76

78

81
82
83

85
86
87
88
89

91

.23268
.22039
.21143
.20258
.19306
.18498
.17679
.16896

.1626]
.15612
.14767
.14094
.13752
.13066
.12485
.11998
.11542
.10902
.10470
.10088
.09539
.09362
.09105
.08638
.08225

.07622
.07025
.06875
.06756
.06454
.06078
.06070
.05545
.05594

.0509]
.05016
.04839
.04579
.0462]
.04498

.0390?
.03775
.03679
.03375
.03275
.03066
.02928

.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100

.57000
.58000
.59000

.61000

.62000

.63000

.65000

.67000

.70000
.71000
.72000
.73000
.74000
.75000
.76000
.77000
.78000

.81000
.82000
.83000
.85000
.87000
.88000

.91000

123

93

95

98

100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138

.03027
.02793
.02685
.02686

.02356
.02230
.02370
.02255
.01957

.0204]

.01740
.01820
.01962
.0169]
.01413
.01500
.01569
.01544
.01405
.0126]
.01318
.01175
.01292
.0102]
.01032
.01019
.00862
.01052
.00970
.00939
.00699
.00638
.00719
.0065]
.00632

.00662

.00624
.00516
.0047]
.00319
.00419
.00503
.0058]

.00100
.00100
.00100
.00100
.00100
.00100
.00100

.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100
.00100

.92000
.93000

.97000
.98000

1.00000
21.01000
1.02000
1.08000
11¥KXXI
L05000
1.06000
1.07000
IJNMXX)
1.09000
1.10000
1.11000
L12000
L13000
1J4000
L15000
L16000
1J7000
1J8000
L19000
L20000
L21000
L22000
L23000
L24000
L25000
L26000
L27000
L28000
L29000
1.30000
1.31000
L32000
1.33000
11.34000
1.35000
1.36000
1.37000

124

APPENDIX E

THE SUBROUTINES MODEL AND SENSE FROM NLINA.FOR
MODIFIED FOR THE DETERMINATION OF THE
OPTIMAL TREATMENT TIME

These subroutines povide the option of using prior information obtained from a
previously performed procedure with the same tumor radius. A second root-finding subroutine
(ZBRENTZ) is included to solve equations (3.70) and (3.74) for the value of 1.... which is a

function of the estimated value of re at each iteration.

SUBROUTINE MODEL
C THIS SUBROUTINE IS FOR CALCULATING ETA. TIE MODEL VALUE
IMPLICIT REAL‘8 (A-H.O-Z)
DIMENSION T(35W,5).Y(35W),SIGZ(35W).B(5).Z(5).
+A(5).BS(5).VINV(5.5).EXTRA(20)
DIMENSION P(5.5).PS(5,5)
C WRITTEN BY JAMES V. BECK
C MODIFED BY LESLE A. SCOTT
C INTHIS PROGRAM,TIETA=ETA ANDETA=BETA
DOUBLE PRECISION KSL. Q. L. P]. ALPHSL. ALPHS
DOUBLE PRECISION ERFC. ZBRENT]
DOUBLE PRECISION ALPHSR. PISR. X. XX. X2. X2ALPH.
+ XALPH. XXALPH. BETA. BETAS. BET A2. ETA, NEWBETA.
+ NEWBE'TS. NEWBETZ. TIME, TIMEC. RADET A1, ETA2

COMMON SIGZ.T.Z.BS.I.ETA.PS.P.B.A.Y.MODL.VII~IV.NP
+.EXTRA
COMMON/MOD/AAJ‘L
COMMON/PROP/KSL. Q. ALPHSL, L. PI. ALPHS
COMMON/TM TIME
COMMON/I‘C/IIMEC
C
EXTERNAL ERFC. ZBRENT]

125

126

KSL = 1.000
ALPHSL = 1.000
L = -1W.000

Q = -1.000

ALPHS = 1.000
RAD = T(I.1)

P1 = DACOS(-1.00)
TIMEC = BS(1)

0

ALPHSR = ALPHSL“0.5D0
PISR = (PI**0.SDO)/2.D0
X IS THE CALCULATED VALUE FOR LAMBDA
ZBRENT is a root finding subroutine used to solve the transcendental equation for the
freezing front location. See Numerical Recipes by Press et al.. Cambridge University
Press. New York. New York, 1986.
X = ZBRENT1(1.0D-4.2.D0.0.001D0)
XX = X*X
X2 = 2.DO*X
X2ALPH = X2*ALPHSR
XALPH = X'ALPHSR
XXALPH = XX*ALPHSL
C TIME = rm
TIME = ZBRENT‘ZCI'IMEC+1.0D-9.TIMEC + 0.010D0.0.001D0)

0000

BETA = RAD/((4.000“'ALPHS*'TI]\E)“0.500)
NEWBETA = RAD/((4.000*ALPHS*(IIME - TIMEC))**0.5000)
C
BETAS = BETA‘BETA
BETA2 = 2.00‘BETA
NEWBETS = NEWBETA‘NEWBET A
NEWBETZ = 2.000‘NEWBETA
C TO INCLUDE PRIOR INFORMATION FROM TIE SAME RADIUS
IF (I .EQ. 1)THEN
ET A] = BS(1)
ELSE
E (BETA.LT.X)TI-IEN
ETA] = 1.000 - Q*(EXP(-BETAS)/BETA2 - EXP(-XX)IX2 -
+ PISR*(ERFC(BETA) - ERFC(X)»
ELSE
E(BETA.GT.X)TIEN
ETA] = (EXP(-ALPHSL*BETAS)/(BETA2*ALPIISR) -PISR*
ERFC(AIJHSR‘BETA))I(EXP(—XXALPH)/X2ALPH

.9.

+ -PISR*ERFC(XALPII))
ELSE
E'TA] = 1.000
ENDE
ENDE

ENDE

127

E (I .EQ. 1) TIEN
ET A2 = 0.000
ELSE
E (NEWBET A .LT. X) TIEN
ETA3 = EXP(-NEWBETS)/NEWBET2
ETA4 = EXP(-XX)/X2
ET A5 = ERFC(NEWBET A)
ET A6 = ERFC(X)
ETA7 =1 - Q‘(ETA3 - ETA4 - PISR*(ETA5 - E'TA6))
ET A2 = ET A7
ELSE
E (NEWBET A .GT. X) TIEN
ETA2 = (EXP(-ALPHSL*NEWBETS)/(NEWBET2*ALPHSR)
+ -PISR"ERFC(ALPHSR“NEWBETA))/
+ (EXP(-XXALPH)/X2ALPH -PISR"ERFC(XALPH))
ELSE
ETA2 = 1.000
ENDE
ENDE
ENDE
ETA = ETAI-ETAZ
RETURN
END
C CALCULATION OF LAMBDA FROM FUNCTION BRENT
DOUBLE PRECISION FUNCTION BRENT1(X1, X2. TOL)

PARAMETER(ITMAX=1W. EPS= 3.0E-8)

DOUBLE PRECISION A. B. C. D. E. FA. FB. FC

DOUBLE PRECISION TOL]. TOL. X1. X2, XM
DOUBLE PRECISION P, Q. R. S. FUNCL]

EXTERNAL FUNCL]
A=X1

=X2
FA=FUNCL1(A)
FB=FUNCL1(B)

IF(FB‘FA .GT. 0.000) PAUSE ’ROOT MUST BE BRACKETED FOR
+ BRENT 1.‘
FfiFB
DO 15 ITER=LITMAX
E(FB"FC .GT. 0.000) TIEN
C=A
FC=FA
D=B-A
E=D
ENDIF
E(ABS(FC) .LT. ABS(FB)) TIEN

15

A=B
B=C
C=A
FA=FB
FB=FC
FC=FA
ENDE

128

TOL1=2.000*EPS*ABS(B)+0.500*TOL

XM=0.500*(C-B)

E(ABS(XM) .LE. TOL] .OR. FB .EQ. 0.000)TIEN

BRENT1=B
RETURN
ENDE

1F(ABS(E) .GE. TOL] .AND. ABS(FA) .GT. ABS(FB)) TIEN

S=FBIFA

IF(A .EQ. C)TIEN
P=2.000*XM*S
$1.000 - S

ELSE
Q=FAIFC
R=FBIFC

P-—S*(2.ODO*XM*Q"'(Q-R) - (B-A)‘(R-1-ODO))
Q=(Q-l.ODO)*(R-1.0DO)"(S-].ODO)

ENDE

E(P .GT. 0) Q = -Q
P=ABS(P)

E(2.000‘P .LT. MTN(3.000‘XM*Q-ABS(TOL1*Q).ABS(E*Q)))TIEN

E=D
0=P/Q
ELSE
D=XM
E=D
ENDE
ELSE
D=XM
B0
ENDE
A=B
FA=FB
E(ABS(D) .GT. TOL1)T'IEN
B=B+D
ELSE
B=B+SIGN(TOLI .XM)
ENDE
FB=FUNCL1(B)
CONTINUE

C

+

129

PAUSE ’BRENT] EXCEEDING MAXIMUM ITERATIONS.’
BRENT 1=B

RETURN

END

DOUBLE PRECISION FUNCTION ERFC(X)

DOUBLE PRECISION A1. A2. A3. A4. A5. P. T. X
A1=0.25482959200
A2=-0.28449673600
A3=1.42141374100
A4=-1.45315202700
A5=1.06140542900
P=0.327591 100
T=1.000/(1 .000-1-P‘X)
ERFC=(A1*T+A2*T“2.000+A3‘T“3.000-0-A4*T"*4.000+A5*T“5.000)
I"EXP(-X"""2.000)
RETURN
END

DOUBLE PRECISION FUNCTION FUNCL](X)

DOUBLE PRECISION KSL. Q. ALPHSL. L. X. P]. ERFC. ALPHS
DOUBLE PRECISION EXPX. EXPXASL. XX2. RATIO
COMMON/PROP/KSL. Q. ALPHSL. L. P]. ALPHS

EXTERNAL ERFC

EXPX=EXP(-X*X)

EXPXASL.-=EXP(-X"X*ALPIISL)

XX2=X*X*2.000

RATIO=EXPXASLI(XX2*ALPHSL“0.500)

FUNCL]=KSL"'Q*EXPX/XX2 - RATIO/(RATIO-(PI”0.500/2)
*ERFC(ALPHSL“‘0.500*X)) -L"'X

RETURN
END

SUBROUTINE SENS

TIIIS SUBROUTINE IS FOR CALCULATING TIE SENSITIVITY COEFFICIENTS

IMPLICIT REAL‘8 (A-H.O-Z)

DIMENSION T(35W.5).Y(35W),SIGZ(35W).B(5).
+Z(5).A(5).BS(5).V1NV(5.5).EXTRA(20)

DIMENSION P(5.5).PS(5.5)

+
4.

DOUBLE PRECISION KSL. Q. ALPHSL. L. P]. ALPHS
DOUBLE PRECISION ERFC. BRENT]. ZBRENTZ
DOUBLE PRECISION ALPHSR. PISR. X. XX. X2. X2ALPH.
XALPH. XXALPH. BET A. BET AS. BET A2. NEWBETA.
NEWBETS. NEWBETZ. DER7. SENSE]. TIME. TIMEC.

n noon

4.

130
DER3. DER6. RAD

COMMON SIGZ.T.Z.BS.I.ETA.PS.P.B.A.Y.MODL.VINV.NP
+.EXTRA

COMMON/MOD/AAJ'L

COMMON/PROP/KSL. Q. ALPHSL. L. P]. ALPHS
COMMON/TC/TIMEC

COMMON/TIMI TIME

COMMON/EQ/X. RAD

+
+
+
+

+

EXTERNAL ERFC. ZBRENT]. ZBRENTZ
KSL = 1.000
ALPHSL = 1.000
L = -100.0D0
Q = -l.0DO
ALPHS = 1.0D0
P1 = DACOS(-1.D0)
RAD = T(I.1)
TIMEC = BS(1)
VARIABLES DECLARED
ALPHSR = ALPHSL**O.SDO
PISR = (PI**O.5D0)/2.DO
X IS THE CALCULATED VALUE FOR LAMBDA
ZBRENTisamotfindingsubroutineusedtosolvedrenanscendentalequafion forthe
freezing from location. See Numerical Recipes by Press et al.. Cambridge University
Press. New York. New York. 1986.
X = ZBRENTI(1.0D-4.2.D0.0.001D0)
WRITE(*.")'X = ‘.X
XX = X‘X
X2 = 2.0D0*X
X2ALPH = X2‘ALPHSR
XALPH = X'ALPHSR
XXALPH = XX*ALPIISL

TIME = ZBRENT2(TIMEC+L00-9.TTMEC + 0.01000.1.00-6)
BETA = RAD/((4.000‘ALPHS‘TIME)“0.500)
NEWBETA = RAD/((4.000*ALPHS"(TIME - TIMEC))**0.5000)
BET AS = BE'TA‘BETA
BET A2 = 2.000‘BET A
NEWBETS = NEWBETA‘NEWBET A
NEWBETZ = 2.000‘NEWBET A

DER3 = Q’EXP(-NEWBETS)/(2.000"NEWBETS)

DER6 8 ((--2.000“NEWBETS’ALPHSL“EXP(-A1..PHSL"l
NEWBETS) - EXP(-ALPHSL"NEWBETS))/
(2.000‘NEWBETS‘ALPHSR) + ALPHSR
*EXP(-ALPHSL*NEWBETS))/((EXP(-ALPHSL*XX))

l(2.000"X*ALPHSR) - PISR’ERFC(ALPHSR‘X))
DER7 = RAD/(4.000*(ALPHS“0.500)*((TIME - TIMEC)"
0.000/2.011))»

+

131

E (1 .EQ. 1) THEN
SENSE] = 1.000
ELSE
E (NEWBET A .LT. X) TIEN
SENSE] = -DER3*DER7
ELSE
E (NEWBETA .GT. X) TIEN
SENSE] = -DER6"'DER7
ELSE
SENSE] = 0.000
ENDE
ENDE
ENDE
Z(1) = SENSE]
RETURN
END

DOUBLE PRECISION FUNCTION BRENT2(X1. X2. TOL)

PARAMETER(ITMAX=1W. EPS= 3.0E-8)

DOUBLE PRECISION A. B. C. D. E. FA. FB. FC

DOUBLE PRECISION TOL]. TOL. X1. X2. XM

DOUBLE PRECISION P. Q. R. S. FUNCLZ. X. TIMEC. RAD

COMMON/TC/TIMEC
COMMON/EQ/X. RAD
EXTERNAL FUNCLZ
A=X1
B=X2
FA=FUNCL2(A)
FB=FUNCL2(B)

E(FB‘FA .GT. 0.000) PAUSE ’ROOT MUST BE BRACKETED FOR
BRENTZ.’
FC=FB
DO 15 TIER=LITMAX .
E(FB‘FC .GT. 0.000) TIEN
C=A
FC=FA
DEB-A
E=D
ENDE
E(ABS(FC) .LT. ABS(FB)) TIEN
A=B
B=C
C=A
FA=FB
FB=FC

15

4..
+

132

FC=FA
ENDIF
TOLl=2.0DO*EPS*ABS(B)o-O.SDO'TOL
XM=O.5DO*(C-B)
E(ABS(XM) .LE.TOL1 .OR. FB .EQ. 0.0DO)THEN
ZBRENT2=B
RETURN
ENDIF
1F(ABS(E) .GE.TOL1 .AND. ABS(FA) .GT. ABS(FB)) THEN
S=FB/FA
IF(A .EQ. C)THEN
P=2.0D0“XM"S
Q=l.ODO - s
ELSE
Q=FAIFC
R=FBIFC
P=S*(2.0DO*XM*Q*(Q-R) - (B-A)*(R-l.OD0))
Q=(Q-r.0D0)*(R-1.0D0)*(s-1.0D0)
ENDIF
IF(P .GT. 0) Q .. -Q
P=ABS(P)
II=(2.0D0*P .LT. MN(3.0D0*XM*Q-ABS(‘TOL1*Q).ABS(E*Q)))THEN
E=D
D=P/Q
ELSE
D=XM
E=D
ENDIF
ELSE
D=XM
E=D
ENDIF
A=B
FA=FB
1F(ABS(D) .GT. TOL1)TI-IEN
B=B+D
ELSE
=B+SIGN(TOL1.XM)
ENDIF
FB=FUNCL2(B)
CONTINUE
PAUSE 'ZBRENTZ EXCEEDING MAXIMUM ITERATIONS.’
ZBRENT2=B
RETURN
END
DOUBLE PRECISION FUNCI'ION FUNCLZCI'IME)
DOUBLE PRECISION KSL. Q. ALPHSL. L. x. P1. ERFC. ALPHS.
DER]. DER2. DER3. DER4. DERs, DER6.
RAD. TIME. BETA. NEWBETA. TIMEC. FUNCLZA. FUNCLZB.

133

+ NEWBETZ. NEWBETS. BET A2. BET AS. PISR. ALPHSR. X2.
+ XX. X2ALPH. XALPH. XXALPH
EXTERNAL ERFC
COMMON/PROP/KSL. Q. ALPHSL. L. PI. ALPHS
COMMON/TC/TIMEC
COMMON/EQ/X. RAD
ALPHSR = ALPHSL”0.500
PISR = (PI“0.500)/2.00
XX = X‘X
X2 = 2.00‘X
X2ALPH = X2‘ALPHSR
XALPII = X‘ALPHSR
XXALPH = XX‘ALPHSL
BETA = RAD/((4.000"ALPHS‘TIME)“0.500)
NEWBET A = RAD/((4.000‘ALPHS‘CTIME - TIMEC))""0.5000)
BET AS = BETA‘BETA
BET A2 = 2.000‘BET A
NEWBETS = NEWBET A‘NEWBET A
NEWBET‘2 = 2.000‘NEWBET A
DER] = Q‘EXP(-BETAS)/(2.000"BETAS)
DER2 = ~RAD/(4.000‘(ALPHS“0.5000)*(TIME**(3.000/2.000)))
DER3 = Q‘EXP(-NEWBETS)I(2.000"NEWBETS)
DER4 = -RAD/(4.000*(ALPHS“0.5000)"((TIME-TIMEC)“
+ (3.000/2.000)))
DER5 = ((-2.000’BETAS‘ALPHSL‘EXP(-ALPIISL*BET AS) -
+ EXP(-ALPHSL"'BETAS))/(2.000"BETAS*ALPHSR)
+ + (ALPHSR*EXH-ALPHSL‘BETAS)))/(EXP(-ALPHSL
+ *XX)/(2.000‘X*ALPHSR) - (PISR‘ERFC(ALPHSR*X)))
DER6 = ((-2.000*NEWBETS*ALPHSL‘EXP(-ALPIISL*NEWBETS) -
+ EXP(-ALPIISL*NEWBETS))/(2.000"NEWBETS*ALPHSR)
+ + (ALPHSR*EXP(-ALPHSL*NEWBETS)))/(EXP(-ALPHSL
+ I"XXI/(2.000"‘X"ALPHSR) - (PISR‘ERFC(ALPHSR"X)))
E (BETA .LT. X) TIEN
FUNCLZA = DERI‘DERZ
ELSE
E (BETA .GT. X) TIEN
FUNCLZA = DER5‘DER2
ENDE
ENDE
E (NEWBET A .LT. X) TIEN
FUNCLZB I: DER3’DER4
ELSE
E (NEWBETA .GT. X) TIEN
FUNCL2B = DER6’DER4
ENDE
ENDE
FUNCL2 = FUNCL2A - FUNCL2B
RETURN
END

APPENDIX F

THE FORTRAN PROGRAM MODC.F OR

This program. MODCFOR. is used to calculate dimensionless temperatures at

corresponding radius locations and times for a given cryosurgical treatment time.

PROGRAM MODC

THIS PROGRAM IS DESIGNED FOR CALCULATING TIE TEMPERATURE
AT A GIVEN R LOCATION AND AT A GIVEN TIME.
WRITTEN BY LESLE A. SCOTT
IMPLICIT REAL‘8 (A-H.O-Z)
DOUBLE PRECISION TIME. TIMEC. R. DELTAT. DELTAR.
+ TIETA. BETA. RINT. TINT
DIMENSION ETA(10). BET(10)

0000

C
COMMON TIETA. BETA
COMMON/EQ/TIME. R. TIMEC
OPEN(UNIT = 14. FILE='TEMPC1.DAT’. ST ATUS=”UNKNOWN")
OPEN(UNIT = 12. F1LE=’BETA.DAT’. ST ATUS="UNKNOWN")

P1 = DACOS(-1.000)
TIMEC = 0.185000
DELTAT = 0.(XX)500
DELTAR = 01an 00
R = 0.1000
RINT = 0.099800
TINT = 0.18000
NR = 5
IR 8 1
NT = 40
WRITE(14.5)(RINT+II"'DELTAR. II = IR. NR-l-TR)
5 FORMAT(9X. 6(1X.F8.4))
DO I = LNT
TIME = TINT+I*DELTAT

134

135

00 II = IR. NR-I-IR
R = RINT+II"DELTAR
CALL MODEL
ET A(II) = TIETA
BET(II) = BETA
C WRITE(12.")BE‘T(II). ET A(II)
ENDDO
WRITE(14.10) TIME. (ET A(TI). II = IR, NR+IR)
WRIIE(14.1 1)(BET (II).II=TR.NR+IR)
10 FORMAT(IXF8.4.6(1X.F8.5))
1 1 FORMAT(7X.6(1X.F8.5))
ENDDO
STOP
END
SUBROUTINE MODEL
C
DOUBLE PRECISION KSL. Q. L. PI. ALPHSL. ALPHS
DOUBLE PRECISION ERFC. BRENT
DOUBLE PRECISION ALPHSR. PISR. X. XX. X2. X2ALPH.
XALPH. XXALPH. BETA. BET AS. BET A2. NEWBET A.
+ NEWBETS. NEWBETZ. ST. TIMEC. TIETA. R. TIME.
+ TIETAZ
COMMON/PROP/KSL. Q. ALPHSL. L. PI. ALPHS
COMMON TIETA. BETA
COMMON/EQ/TIME. R. TIMEC

+

C
EXTERNAL ERFC. ZBRENT
KSL = 1.0D0
ALPHSL = 1.000
L = -lO0.0DO
Q = -l.0D0
ALPHS = 1.000
P1 = DACOS(-I.D0)
ALPHSR = ALPHSL**0.5D0
PISR = (PI**0.5D0)/2.D0
X 18 THE CALCULATED VALUE FOR LAMBDA
ZBRENTisamotfindmgsubmufimusedmsolvemenanscenduualequafionforme
freezing from location. See Numerical Recipes by Press et al.. Cambridge University
Press. New York. New York. 1986.
X a ZBRENT(1.0D-4.2.D0.0.00IDO)
XX 2 X*X
X2 = 2.D0‘X
X2ALPH = XZ‘ALPHSR
XALPII s X'ALPHSR
XXALPH a XX*ALPHSL

0000

BETA = R/((4.000*ALPHS"TIME)“0.500)
ST = X*2.000"‘((ALPHS"TII\E)“0.5000)
BET AS = BET A‘BET A

C

+

136

BET A2 8 2.TXI"BETA
E (BETA .LT. X) TIEN
TIETA 8 ] - Q‘(EXP(-BET AS)/BET A2 - EXP(-XX)/X2 -
PISR*(ERFC(BET A) - ERFC(X)»
ELSE
E (BETA .GT. X) TIEN
TIETA = (EXP(-ALPHSL"BETAS)/(BETA2*ALPIISR) -PISR
‘ERFC(ALPHSR"BETA))/(EXP(-XXALPH)/X2ALPH
~PISR‘ERFC(XALPH))
ELSE
TIETA = 1.000
ENDE
ENDE
E (TIME .GT. TIMEC) TIEN
NEWBETA = R/((4.000"ALPHS"(TIME - TIMEC»“0.5000)
NEWBETS = NEWBETA‘NEWBETA
NEWBETZ = 2.000‘NEWBETA
E (NEWBETA .LT. X) TIEN
TIETA2 = 1 - Q*(EXP(-NEWBETS)/NEWBET2 - EXP(-XX)/X2
- PISR*(ERFC(NEWBETA) - ERFC(X)»
WRITE(“.‘)TIETA2
ELSE

E (NEWBETA .GT. X) TIEN
TIETA2 = (EXP(-ALPHSL*NEWBETS)/(NEWBET2*ALPHSR)
-PISR*ERFC(ALPHSR*NEWBETA))I
(EXP(-XXALPH)/X2ALPH -PISR"ERFC(XALPII))
ELSE
TIETA2 = 1.000
ENDE
ENDE
TIETA = TIETA - TIETA2
ENDE
RETURN
END
CALCULATION OF LAMBDA FROM FUNCTION BRENT
DOUBLE PRECISION FUNCTION BRENT(X1. X2. TOL)
PARAMETER(ITMAX=1W. EPS= 3.0E-8)
DOUBLE PRECISION A. B. C. D. 5. FA. FB. FC
DOUBLE PRECISION TOL]. TOL. X1. X2. XM
DOUBLE PRECISION P. Q. R. S. FUNCL
EXTERNAL FUNCL
A=X1
B=X2
FA=FUNCL(A)
FB=FUNCL(B)

E(FB‘FA .GT. 0.000) PAUSE 'ROOT MUST BE BRACKETED FOR
BRENT.‘
FC=FB

137

00 15 ITER=1 .ITMAX
E(FB‘FC .GT. 0.000) TIEN
C=A
FC=FA
D=B-A
E=D
ENDE
E(ABS(FC) .LT. ABS(FB)) TIEN
A=B
B=C
@A
FA=FB
FB=FC
FC=FA
ENDE

TOL1=2.000“EPS*ABS(B)+0.500"TOL
XM=0.500‘(C-B)
E(ABS(XM) .LE. TOL] .OR. FB .EQ. 0.000)TIIEN
BRENT=B
RETURN
ENDE

E(ABS(E) .GE. TOL] .AND. ABS(FA) .GT. ABS(FB)) TIEN
S=FBIFA
E(A .EQ. C)TIEN
P=2.000"XM*S
Q=1.000 - S
ELSE
Q=FAIFC
R=FBIFC
P=S*(2.0N*XM*Q‘(Q-R) - (B-A)’(R-1.000))
Q=(Q- 1.000)‘(R-1.000)"(S-1.000)
ENDIF
E(P .GT. 0) Q = -Q
P=ABS(P)
E(2.000"‘P .LT. W(3.0m‘)CM*QABSCTOL1*Q)ABS(E*Q)))TIEN
E=D
D=P/Q
ELSE
D=XM
E=D
ENDE

ELSE

D=XM
E=D

ENDE
A=B
FA=FB

138

E(ABS(D) .GT. TOL1)TIEN
B=B+D
ELSE
B=B+SIGN (T OLl .XM)
ENDE
FB=FUNCL(B)
CONTINUE

PAUSE ‘BRENT EXCEEDING MAXIMUM ITERATIONS.’
BRENT =B

RETURN

END

DOUBLE PRECISION FUNCTION ERFC(X)

DOUBLE PRECISION A]. A2. A3, A4. A5. P. T. X
A1=0.25482959200
A2=-0.28449673600
A3=1.42141374IIX)
A4=-1.45315202700
A5=L06140542900
P-0.327591100
T=1.000/(1 .000+P"'X)
ERFC=(A1*T+A2"T**2.0m+A3"'T“3.000+A4*T“4.000+A5*T“5.000)
+ ‘EXP(-X‘”"2.000)
RETURN
END

DOUBLE PRECISION FUNCTION FUNCL(X)

DOUBLE PRECISION KSL. Q. ALPHSL. L. X. PI. ERFC. ALPHS
DOUBLE PRECISION EXPX. EXPXASL. XX2. RATIO
COMMON/PROP/KSL. Q. ALPHSL. L. PI. ALPHS

EXTERNAL ERFC
EXPX=EXP(-X"X)
EXPXASL=EXP(-X"'X*ALPHSL)
XX2=X*X*2.0DO
RATTO=EXPXASU(XX2‘ALPHSL”0.500)
FUNClzKSL‘Q‘EXPX/XXZ - RATIO/(RATIO-(PI“0.500/Z)
+ 1“ERFC(ALPIISL"“0.500“'X)) -L"X
RETURN
END

139
ThisfilerepresentstheoutputfilefromthepmgramMODCFOR. Thefirstrowisthe

radiusvalues.thefirstentryistbesecondrowisthetime. Theremainingentriesinthesecond

rowaredimensionlesstemperatures. Thisaltemafingpattemisrepeatedtiuougboutthelist.

.1805
.1810
.1815
.1820
.1825
.1830
.1835
.1840
.1845
.1850
.1855
.1860
.1865
.1870
.1875
.1880
.1885
.1890
.1895

.19W

.0999
2.33546
.11757
2.34126
.1 1741
2.34706
.1 1725
2.35285
.1 1708
2.35863
.1 1692
2.36441
.11676
2.37017
.1 1661
2.37593
.1 1645
2.38168
.1 1629
2.38743
.11613
2.39311
.11597
2.39721
.1 1582
2.39837
.11566
2.39724
.1 1551
2.39464
.11535
2.391 13
.1 1520
2.38709
.1 1505
2.38272
.1 1490
2.37819
.1 1474
2.37359
.11459

.1(XX)
2.33126
.1 1769
2.33706
.1 1753
2.34285
.11736
2.34864
.1 1720
2.35441
.1 1704
2.36018
.1 1688
2.36594
.1 1672
2.37170
.1 1656
2.37744
.1 1641
2.38318
.1 1625
2.38885
.11609
2.39296
.11593
2.39414
.11578
2.39303
.1 1562
2.39045
.11547
2.38696
.1 1532
2.38294
.1 1516
2.37859
.1 1501
2.37408
.1 1486
2.36949
.1147]

.1W1
2.32708
.1178]
2.33287
.11764
2.33866
.11748
2.34443
.11732
2.35020
.11716
2.35597
.117W
2.36172
.11684
2.36747
.11668
2.37321
.11652
2.37894
.11636
2.38461
.1162]
2.38872
.11605
2.38991
.11589
2.38882
.11574
2.38627
.11559
2.38280
.11543
2.37879
.11528
2.37447
.11513
2.36997
.1149?
2.36539
.11482

.1W2
2.32290
.11792
2.32869
.11776
2.33447
.11760
2.34024
.11744
2.346W
.11728
2.35176
.1171]
2.35751
.11696
2.36325
.11680
2.36898
.11664
2.37471
.11648
2.38037
.11632
2.38449
.11617
2.38570
.1160]
2.38463
.11586
2.38209
.11570
2.37865
.11555
2.37466
.11539
2.37035
.11524
2.36587
.11509
2.36131
.11494

.1W3
2.31873

.11804

2.3245 1

.11788

2.33029

.11772

2.33605

.11755

2.34181

.11739

2.34756

.11723

2.35330

.11707

2.35904

.1169]

2.36477

.11675

2.37049

.11660

2.37615

.11644

2.38027

.11628

2.38150

.11613

2.38045

.11597

2.37793

.11582

2.3745 1

.11566

2.37053

.1155]

2.36624

.11536

2.36178

.11520

2.35723

.11505

. 1W4
2.31457
.11816
2.32035
.118W
2.32611
.11783
2.33187
.11767
2.33763
.1175]
2.34337
.11735
2.34911
.11719
2.35484
.11703
2.36056
.11687
2.36628
.1167]
2.37193
.11656
2.37606
.11640
2.37730
.11624
2.37627
.11609
2.37377
.11593
2.37037
.11578
2.36642
.11562
2.36214
.11547
2.35769
.11532
2.35316
.11517

APPENDIX G

THE FORTRAN PROGRAM RAD.FOR

This program. RAD.FOR. was written to read a file Of data. TEMPDAT. and to add

random errors to the fourth column of that data file. This simulates random measurement errors

present in the radius measurement data. Both TEMP.DAT and the output file. TEMPCDAT. are

used as input to the program NLINA.FOR for the determination of the Optimal treatment time. t,

C

PROGRAM RAD

C WRITTEN BY DEBBE MONCMAN

COMMON/RAND/NT. STDDV
COMMON NDAT.NPTS.C3
DOUBLE PRECISION B(4).X(5)
DIMENSION DATA(2(XXX))

INTEGER C2.NT.C5.C6.C7.I.C3.C4. C8

OPEN(UNTT=14. FILE = ‘TEMP.DAT’. STATUS = ’UNKNOWN’)
READ(14."') NT.C2.C3.C4.C5.C6

READ(14.*) (B(I).I=1.C2)

OPEN(UNIT815. FILE=°TEMPC.DAT'.STATUS=’UNKNOWN’)
WRITE(ISJ) NT. C2. C3. C4. C5. C6

FORMAT(IX. 6(15.1X))

WRITE(15.6) (B(I).I=1.C2)

FORMAT(lX.4(E14.8.1X))

CALL RANDOM (DATA)

DO 10. I = 1. NT
READ(14.*)II. TEMP. SIGMA. “(Db-"LC”
WRITE(15.7)II.TEMP.SIGMA.(X(I)-I-DATA((I-1)‘C3+J).J=1.C3)
FORMAT(IXJ5.1X.5(E14.8.1X))

CONTINUE

READ(14.")C7. C8

140

000 0 0 0099 0

000

11
12

141

WRITE( 15.8) C7. C8
FORMAT( 1X.I5)
CLOSE( 14)

CLOSE( 15)

STOP

END

SUBROUTINE RANDOM (DATA)

INTEGER C3

COMMON/RAND/NT. STDDV

COMMON NDATNPT‘S.C3

PARAMETER(PI=3. 14159265 .NP'I‘Sa4.NBIN-=1000.NDAT=NP'TS+NBIN)

PARAMETER(PI=3. 14159265.NBIN=10W)

SEE Numerical Recipes by Press. Flamtery. Teukolsky and Vetteriing.
Cambridge Press. 1986 about page 192

Modified by JV Beck. Michigan State University.

E-mail address: 22427jvb@ibm.cl.msu.edu

DIMENSION DATA(2(I)00)

CHARACTER‘SO FOUT
NPTS = NT*C3
IDUM IS SEED. SET TO ANY NEGATIVE NUMBER TO INITIALIZE OR
REINITIALIZE.
IDUM=~5

WRITE(*.*)’ Enter the number of points ’

READ(*.'*)NPTS

WRITE(*.*)’ENTER TIE SEED NUMBER (-)’

READ(*.*)IDUM

NDAT=NPTS+NBIN

WRITE(*.*)’ GIVE TIE STANDARD DEVIATION’

READ(*.")STDDV

WRITE(*.*)’Give the name of the output file’
READ(*.’(A80)’)FOUT
OPEN(13. FILE=FOUT)

RHON=0.0

RHOD=0.0

WRITE(*.*)’ I RAND. NO.’

DO 500 IDUMI=1.1

DATA(1)=GASDEV(IDUM)*STDDV

WRITE(*.100)1.DATA(1)

WRITE(13.1W)1.DATA(1)

DO 11 1=2.NP'TS
DATA(I)=GASDEV(IDUM)*STDDV
RHON=RHON+DATA(1-l)*DATA(1)
RHOD=RHOD+DATA(I)"DATA(I)

WRITE(*.100)I.DATA(I)

WRITE(13.1W)I.DATA(I)
CONTINUE
CONTINUE
RHO=RIIONIRHOD

142

CCC WRITE(*.‘(1X.AI)’) ’Descriptors of a gaussian distribution’

CALL MOMENT (DATAJ-IAVEADEVSDEVNARRHO)

5W CONTINUE

C

WRIT’E(*.’(1X.’I‘29.A.T42.A/)’) ’ Values of qrantities’.’ ’
WRITE(‘J'Y Values of quantities’
WRITE(*.’(1X.T29.A.T42.A/)’) ’ Sample ’JExpected’
WRITE(*.’(1X.A.T25.2F12.4)’) ’Mean :’.AVE.0.0
WRITE(*.'(1X.A.T25.2F12.4)’) ’Average Deviation :’.ADEV.STDDV
WRITE(*.’(1X.A.T25.2F12.4)’) ’Standard Deviation :’.SDEV.STDDV
VARTH=STDDV*STDDV
WRITE(‘.’(1X.A.T25.2F12.4)’) ’Variance :’.VAR.VARTH
WRI’I'E(*,’(1X.A.T25.F12.4)’)’Est. Correlation Coef.’.RHO
WRITE(*.*)’Aver-age deviation comes from use of absolute values’

1W FORMAT(I]0.F 10.6)

11

12

END
SUBROUTINE MOMENT(DATA.N.AVE.ADEV.SDEV.VAR.RHO)
DIMENSION DATA(2(XXX))
E(N.LE.1)PAUSE ’N must be at least 2’
S=0.
80:0.
SN=0.
DO 11 J=LN
S=S+DATA(J)
E(I .EQ. 1)GOTO 11
SN=SN+DATA(J)*DATA(J-1)
SD=SD+DATA(I)*DATA(I)
CONTINUE
AVE=SIN
ADEV=0.
VAR=0.
DO 12 J=I.N
S=DATA(.I)—AVE
ADEV=ADEV+ABS(S)
P=S’S
VAR=VAR+P
CONTINUE
ADEV=ADEVIN
VAR=VAR/(N-1)
SDEV=SQRT(VAR)
RHO=SN/SD
WRITE(‘J‘YSN SD RIIO’.SN.SD.RHO
RETURN
END
FUNCTION RAN1(IDUM)
DIMENSION R(97)

RETURNS UNEORMLY DISTRIBUTED NUMBERS BETWEEN 0 AND 1

PARANETER (M1=2592W.IA1=7141.IC1=54773.RM1=3.8580247E—6)
PARAMETER (M2=134456.1A2=8121.IC2=28411.RM2=7.4373773E-6)
PARAMETER (M3=243(XX).IA3=4561.IC3=51349)

143

DATA EF MI
E (IDUM.LT.0.0R.EF.EQ.0) TIEN

EF=1

IX1=MOD(IC1-IDUM.M1)

IX1=MOD(IA1*IX1+IC1.M1)

IX2=MOD(IX1.M2)

IX1=MOD(IA1"IX1+IC1.M1)

IX3=MOD(IX1.M3)

DO 11 J=1.97
IX1=MOD(IA1"IX1+IC1.M1)
IX2=MOD(IA2‘IX2+ICZ.M2)
R(J)=(FLOAT(IX1)+FLOAT(D(2)*RM2)"RM1

1 1 CONTINUE
IDUM=1
ENDE
IX1=MOD(IA1"IX1+IC1.M1)
IX2=MOD(IA2‘IX2+IC2.M2)
IX3=MOD(IA3’IX3+IC3.M3)
J=1+(97*IX3)/M3
E(J.GT.97.0RJ.LT. ])PAUSE
RAN 1=R(J)
R(WAT(IX1)+FLOAT(IX2)*RM2)"RM1
C WRl'I'EC“.“)’J.R(I).l?~A1*I1’.J.l?»(J).I?~ANl
RETURN

END
FUNCTION GASDEV(IDUM)
C USES BOX-MULLER TRANSFORMATION FROM UNEORM DISTRIBUTION TO
C NORMAL DISTRIBUTION WITI'I UNIT STANDARD DEVIATION
DATA ISET/O/
E (ISET.EQ.0) TIEN
1 V1=2."RAN1(IDUM)—1.
V2=2."‘RAN1(IDUM)-1.
R=V1**2+V2“2
E(R.GE.1..OR.R.EQ.0.)GO TO 1
FAC=SQRT(-2.*LOG(R)/R)
GSET=V1"'FAC
GASDEV=V2*FAC
ISET=1
ELSE
GASDEVaGSET
ISET =0
ENDE
C WRITE(‘.‘)’IDUM.GASDEV’.IDUM.GASDEV
RETURN
END

as input for NLINA.FOR when exact radius measurement data was used. The first row of
numbersrepresentthenumberofdatapoints,thenumberofparameterstobeestimated.the
number Of independent variables. the maximum munber of iterations to be performed. the model
number.andtheusualprintoutsrespectively. Thesecondrowrepresentstheinitialguessofthe
Optimalneannerucoolirrgtimetobedetermined. T'hefirstcolumnistheindex.thesecond
columnisundesireddimensiomesstempemmres.memirdismestandarddeviafionoffire

measurement errors. and the fourth is the measured values for the radii. without random enors.

addedtotheradiusmeasurementsinthefourthcolumn. ItisalsousedasinputforNLINAFOR.

10

This file. TEMPDAT. represents the input to the program RAD.FOR. It was also used

10.1.1.2W.1.1

0.125000
1 2.39414
2 2.39414
3 2.39414
4 2.39414
5 2.39414
6 0.99467
7 0.99467
8 0.99467
9 0.99467
10 0.99467

0

This file. TEMPC.DAT. is an output file from RADFOR with random measurement enors

1 1 2W

.125(XXXX)E+W

Somqaus-mto—

.239414WE+01
.239414WE+01
.239414WE-l-01
.239414WE+01
.239414WE+01
.99467WOE+W
.99467WOE+W
.99467WOE+W
.99467WOE+W
.99467(XX)E+W

0.W10
0.W10
0.W10
0.W10
0.W10
0.W10
0.W10
0.W10
0.W10
0.W10

1 1

JUNE-02
AWE-02
.11XXX3WOE-02
AWE-02
.1000000015-02
.IW-W
.ImE-OZ
.100000005-02
.1mWE-02
.lmE-OZ

144

0.1000
0.1111)
0.1W0
0.1W0
0.1W0
0.15W
0.15W
0.15W
0.15W
0.15W

.10147277E+W
.99968570E-0]

.999971085-01

.1W73597E+W
.99912973E-01

.15010918E4-W
.14650322E+W
.1490] 197E+W
.1502657654-00
.14862881E+W

APPENDIX H

THE SUBROUTINES MODEL AND SENSE FROM NLINA.FOR
MODIFIED FOR THE DETERMINATION OF THE OPTIMAL
TREATMENT TIME WITH PRIOR INFORMATION
FROM A DIFFERENT RADIUS

Inthemainprogram.theterrnEXTRA(l)issetequalto 1 whenpriorinformationfrom
asingle different radius isused. andsetequal t02 whentwo differentradiiare used. Calculations
performed in the subroutine MODEL are as follows: beginning with the estimate of the treatment
time. I... and the corresponding radius. Rad,. the time that the minimum temperature is achieved.
1...... is calculated. Using these values, the minimum temperature. T... is determined. Using 1'...
and the radius at which the prior information was obtained. Rad, a calculated value of rd is then
Obtained. ThisvalueisusedwiththeactualvalueoftwobtainedfromtheinputfileJnthe
modified sum of squares function given by equation (3.4).

Inthe SENSE mbroutine.thefinitedifferencemethodis usedtodetermine thesensitivity

coefficients when EXTRA(1) is equal to 1 or 2.

SUBROUTINE MODEL
C TIIIS SUBROUTINE IS FOR CALCULATING ETA. TIE MODEL VALUE
IMPLICIT REAL‘8 (A-H.O«Z)
DIMENSION T( 35W.5).Y(35W).SIG2(35W).B(5)1(5).
+A(5).BS(5).VINV(5.5).EXTRA(20)
DIMENSION P(5.5).PS(5.5)
C WRITTEN BY JAMES V. BECK
C MODEED BY LESLE A. SCOTT
DOUBLE PRECISION KSL. Q. L. PI. ALPHSL. ALPHS
DOUBLE PRECISION ERFC. BRENT 1

145

I46

DOUBLE PRECISION ALPHSR. PISR. X. XX. X2. X2ALPH.

+ XALPH. XXALPH. BETA. BETAS. BETA2. ETA. NEWBETA.
+ NEWBETS. NEWBETZ, TIME. TIMEC. RAD. ETAl. ETA2.
+ RAD2

COMMON SIGZ.T.Z.BS,I.ETA.PS.P.B.A.Y.MODL.VINV.NP
+.EXTRA

COMMON/MOD/AAII'L

COMMON/PROP/KSL. Q. ALPHSL. L. PI. ALPHS

C
EXTERNAL ERFC. ZBRENTl

KSL = 1.000

ALPHSL = 1.0D0
L = -1(X).OD0
Q = ~1.0D0

ALPHS = 1.0D0
RAD = T(I.2)
RAD2 = T(I.1)

PI = DACOS(-1.D0)
TIMEC = BS(1)

D

ALPHSR = ALPHSL**O.5D0

PISR = (Pl“0.5D0)/2.D0
X IS THE CALCULATED VALUE FOR LAMBDA
ZBRENTisamotfindingsubmutineusedtosolvetheuamcendentalequationforthe
freezing from location. See Numerical Recipes by Press et al.. Cambridge University
Press. New York. New York. 1986.
X = ZBRENTl(l.OD-4.2.D0.0.001D0)

XX = X*X

X2 = 2.D0"'X

X2ALPH = X2*ALPHSR

XALPH = X'ALPHSR

XXALPH = XX*ALPHSL

0000

TIME = ZBRENTZCI'MEC+1.0D-9.TIMEC + 0.0lOD0.l.OD-6. RAD. TIMEC)

BETA = RAD/((4.0D0‘ALPHS‘TIME)“O.SDO)
NEWBETA = RAD/((4.0D0*ALPHS*(TIME - T'IMEC))“0.SODO)

BET AS = BETA’BETA
BETA2 = 2.D0*BET A
NEWBETS = NEWBET A‘NEWBET A
NEWBETZ = 2.0D0‘NEWBET A
IF (1 .LE. EXTRA(1))THEN
CALL MODEL2(BS(1). ETAl)
ELSE
IF (BET A.LT.X)THEN
ETA] = 1.000 - Q‘(EXP(-BETAS)/BETA2 - EXP(-XX)/X2 -
+ PISR“(ERPC(BET A) - ERFC(X)»

147

ELSE
IF(BE'I‘A.GT.X)THEN
ETA] = (EXP(-ALPHSL*BETAS)I(BETA2"ALPHSR) -PISR"
+ ERFC(ALPHSR‘BETA))/(EXP(-XXALPH)/X2ALPH
+ -PISR*ERFC(XALPH))
ELSE
ETA] = 1.000
ENDIF
ENDIF
ENDIF

IF (1 .LE. EXTRA(1)) THEN
ET A2 = 0.0D0
ELSE
IF (NEWBETA .LT. X) THEN
ETA3 = EXH-NEWBETSVNEWBETZ
ETA4 = EXP(-XX)/X2
ETA5 = ERFC(NEWBET A)
ETA6 = ERFC(X)
ETA7 = l - Q‘CET A3 - ETA4 - PISR“(ETA5 - ETA6))
ET A2 = ET A7
ELSE
IF (NEWBET A .GT. X) THEN
ET A2 = (EXP(-ALPHSL"NEWBETS)/(NEWBET2*ALPHSR)
+ -PISR*ERFC(ALPHSR*NEWBETA))/
+ (EXP(-XXALPH)/X2ALPH -PISR*ERFC(XALPH))
ELSE
ET A2 = 1.0D0
ENDIF
ENDIF
ENDIF
ETA = ETAl-ETA2
IF (I .EQ. 1) THEN
TMIN = ETA
EN DIF
RETURN
END
C CALCULATION OF LAMBDA FROM FUNCTION ZBRENT
DOUBLE PRECISION FUNCTION ZBRENT1(X]. X2. TOL)

000

C
PARAMETER(ITMAX=1W. EPS= 3.0E-8)
DOUBLE PRECISION A. B. C. D. E. FA. FB. FC
DOUBLE PRECISION TOL]. TOL. X1. X2. XM

DOUBLE PRECISION P. Q. R. S. FUNCL]

C
EXTERNAL FUNCL]
A=Xl
B=X2

FA=FUNCL1(A)

148

FB=FUNCLI(B)
C
IF(FB‘FA .GT. 0.0D0) PAUSE ’ROO'I‘ MUST BE BRACKETED FOR
+ ZBRENT] .‘
FaFB
DO 15 ITER=1 .l'l'MAX
IF(FB‘FC .GT. 0.0D0) THEN
C=A
FC=FA
D=B-A
E=D
ENDIF
E(ABS(FC) .LT. ABS(FB)) THEN
A=B
B=C
C=A
FA:
FB=FC
FC=FA
ENDIF

TOL]=2.0DO"EPS*ABS(B)+O.5D0‘TOL
XM=O.5DO*(C-B)
E(ABS(XM) .LE. TOL] .OR. FB .EQ. 0.0D0)TIEN
ZBRENT1=B
RETURN
ENDIF

1F(ABS(E) .GE. TOLl .AND. ABS(FA) .GT. ABS(FB)) THEN
S=FBIFA
IF(A .EQ. C)THEN
P=2.0D0‘XM*S
Q=l.0DO - S
ELSE
Q‘FA/FC
R=FB/FC
P=S"'(2.0DO*XM"Q*(Q-R) - (B-A)*(R-l.0D0))
Q==(Q-].0IX))*(R-].ODO)“(S-l.0m)
ENDIF
IF(P .GT. 0) Q = —Q
P=ABS(P)
IF(2.0D0"P .LT. W(3.0D0"m*Q«ABS(TOLl*Q),ABS(E*Q)))TT-IEN
E=D
D=P/Q
ELSE
D=XM
E=D
ENDIF
ELSE

149

D=XM
E=D

ENDIF

A=B

FA=FB

1F(ABS(D) .GT. TOL])THEN
B=B+D

ELSE
B=B+SIGN(TOL1.XM)

ENDIF

FB=FUNCL](B)

]5 CONTINUE

PAUSE 'ZBRENT] EXCEEDING MAXIMUM ITERATIONS.‘
ZBRENT1=B

RETURN

END

DOUBLE PRECISION FUNCTION ERFC(X)

DOUBLE PRECISION A]. A2. A3. A4. A5. P. T. X
A1=0.254829592D0
A2=-O.284496736D0
A3=l.42]4l 3741D0
A4=-].453152027D0
A5=l .061405429D0
P=0.327591 IDO
T=].0D0/(1.0DO+P"X)
ERFC:(A1*T-r-AZ‘T"*2.0D0+A3"T“3.0D0+A4*T“4.0D0+A5“T“SDDO)
+ *EXP(-X“2.0D0)
RETURN
END

DOUBLE PRECISION FUNCTION FUNCL](X)

DOUBLE PRECISION KSL. Q. ALPHSL. L. X. PI. ERFC. ALPHS
DOUBLE PRECISION EXPX. EXPXASL. XX2. RATIO
COMMON/PROP/KSL. Q. ALPHSL. L. PI. ALPHS

EXTERNALERFC

EXPX=EXP(-X*X)
EXPXASLzEXH-X“X*ALPHSL)
XX2=X"X*2.0D0
RATIO=EXPXASU(XX2*ALPHSL“O.5DO)
FUNCL]=KSL*Q*EXPX/X.X2 - RATTO/(RATTO-(PI**0.5D0/2)
+ *ERFC(ALPHSL“O.5DO*X)) -L"X
RETURN
END

C

150

SUBROUTINE SENS

C THIS SUBROUTTNE IS FOR CALCULATING THE SENSITIVITY COEFFICIENTS

C

O 0000

IMPLICIT REAL‘S (A-H.O-Z)
‘ DIMENSION T(3500.5).Y(3500),8162(3500).B(5),
+Z(5).A(5).BS(5).VINV(5.5).EX'I'RA(20)
DIMENSION P(5.5).PS(5.5)

DOUBLE PRECISION KSL. Q. ALPHSL. L. PI. ALPHS
DOUBLE PRECISION ERFC. ZBRENT]. ZBRENT?
DOUBLE PRECISION ALPHSR. PISR. X. XX. X2. X2ALPH.
XALPH. XXALPH, BETA. BET AS. BETA2. NEWBETA.
NEWBETS. NEWBETZ. DER7. SENSE]. TIME. TIMEC.
DER3. DER6. RAD. RAD2. TMIN. ETA. ET AB

COMMON 8162.11.38.I.E'l‘A.PS.P.B.A.Y.MODL.VINV.NP
+.EX'I'RA
COMMONMOD/AAJ'L

COMMON/PROP/KSL. Q. ALPHSL. L. PI. ALPHS
COMMON/EQ/X
COMMON/RD/RAD. RAD2
COMMON/TEMP/TMIN

EXTERNAL ERFC. ZBRENT]. ZBRENTZ
KSL = 1.0D0
ALPHSL = mm
L = -100.0D0
Q = -l.0D0
ALPHS = 1.000
PI = DACOS(-1.D0)
RAD = T(I.2)
RAD2 = T(I.1)
TIMEC = BS(1)
VARIABLES DECLARED
ALPHSR = ALPHSL“0.5D0
PISR = (Pl“0.5D0)/2.D0
X IS THE CALCULATED VALUE FOR LAMBDA
ZBRENTisarootfindingsubroutineusedtosolvethetranscendentalequationforthe
freezing from location See Numerical Recipes by Press et al.. Cambridge University
Press. New York. New York. 1986.
X = ZBRENT1(1.0D-4.2.D0.0.001DO)
WRIT'E(*.‘)’X = ’.X
XX = X‘X
X2 = 2.0D0*X
X2ALPH = X2*ALPHSR
XALPI-l = X‘ALPHSR
XXALPH = XX'ALPI-ISL

TIME = ZBRENT‘ZCTIMEC+1.0D-9.TTMEC + 0.0]0D0.l.0D-6. RAD. TIMEC)
BETA s RAD/((4.0DO*ALPHS*TIME)“O.5DO)

15]

NEWBET A = RAD/((4.0DO‘AIJ’HS‘CTTME - TIMEC))“0.50D0)
BETAS = BETA‘BETA
BET A2 = 2.0WBETA
NEWBETS = NEWBETA'NEWBET‘A
NEWBET‘Z = 2.0DO"NEWBET A

DER3 = Q’EXP(-NEWBETS)/(2.0D0*NEWBETS)

DER6 = ((-2.0D0“NEWBETS"AI.PHSL"EXP(-ALPHSL"I
NEWBETS) - EXP(-ALPHSL*NEWBE‘TS))/
(2.0D0*NEWBETS*ALPHSR) + ALPHSR
‘EXP(-ALPHSL*NEWBETS))/((EXP(~ALPHSL*XX))

l(2.0DO*X*ALPHSR) - PISR‘ERFC(ALPHSR*X))
DER7 = RAD/(4.0DO"(ALPI-IS“O.5DO)*((TTME - TIMEC)"
+ (3.0D0f2.0D0)))

++++

C
IF (I .LE. EXTRA(])) THEN
CALL MODEL2(BS(1). ETA)
ET AB = ETA
CALL MODEL2(BS(1)*(].OD0+1.0D—]2).ET A)
SENSE] = (ETA - ETAB)/1.0D-12
ETA = ETAB
ELSE
IF (NEWBETA .LT. X) THEN
SENSE] = -DER3"'DER7
ELSE
IF (NEWBET A .GT. X) THEN
SENSE] = -DERG*DER7
ELSE
SENSE] = 0.0D0
ENDIF
ENDIF
ENDIF
20) = SENSE]
RETURN
END

DOUBLE PRECISION FUNCTION ZBRENT2(X]. X2. TOL. RAD. TIMEC)

PARAMETER(ITMAX=1W. EPS= 3.0E-8)

DOUBLE PRECISION A. B. C. D. E. FA. FB. FC

DOUBLE PRECISION TOLl, TOL. X1. X2. XM

DOUBLE PRECISION P. Q. R. S. FUNCLZ. X. TIMEC. RAD

COMMON/EQ/X
EXTERNAL FUNCLZ
A=Xl
B=X2
FA=FUNCL2(A. RAD. TIMEC)
FB=FUNC12(B. RAD. TIMEC)

C

+

152

IF(FB‘FA .GT. 0.0D0) PAUSE ’ROOT MUST BE BRACKETED FOR
ZBRENTZ.’
FC=FB
DO 15 ITER=1JTMAX
IF(FB‘FC .GT. 0.0D0) THEN
C=A
FO-FA
D=B-A
E=D
ENDIF
E(ABS(FC) .LT. ABS(FB)) THEN
A=B
B=C
C=A
FA=FB
FB=FC
FC=FA
ENDIF

TOL1=2.0D0"EPS*ABS(B)+O.5D0*TOL
XM=O.5D0"(C-B)
1F(ABS(XM) .LE. TOL] .OR. FB .EQ. 0.0D0)TIEN
ZBRENTZ=B
RETURN
ENDIF

1F(ABS(E) .GE. TOL] .AND. ABS(FA) .GT. ABS(FB)) THEN
S=FBIFA
IF(A .EQ. C)TI'IEN
P=2.0D0"XM"S
Q=1.0DO - S
ELSE
Q=FAIFC
RSFB/FC
P=S“(2.0W*XM*Q*(Q-R) - (B-A)*(R-I.OD0))
Q=(Q-I.0IX))’(R-l.0m)‘(S-l.0m)
ENDIF
IF(P .GT. 0) Q = -Q
hABSCP)
IF(2.0D0"P .LT. MIN(3.0WXM*Q-ABS(TOL]*Q).ABS(E*Q)))TT'IEN
E=D
DIP/Q
ELSE
D=XM
E=D
ENDIF
ELSE
D=XM

153

E=D
ENDIF
A=B
FA=FB
IF(ABS(D) .GT. TOL])THEN
B=B+D
ELSE
B=B+SIGNCTOL1.XM)
ENDIF
FB=FUNQ.2(B. RAD. TIMEC)
15 CONTINUE
C
PAUSE 'ZBRENTZ EXCEEDING MAXIMUM ITERATIONS.’
ZBRENT2=B
RETURN
END

DOUBLE PRECISION FUNCTTON FUNCLZCT'IME. RAD. TIMEC)

DOUBLE PRECISION KSL. Q. ALPHSL. L. X. PI. ERFC. ALPHS.
DER]. DER2. DERB. DER4. DER5. DERG.
RAD. TIME. BETA. NEWBETA. TIMEC. FUNCLZA, FUNCLZB.
NEWBETZ. NEWBETS. BETA2. BET AS. PISR. ALPHSR. X2.
XX. X2ALPH. XALPH. XXALPH

++++

EXTERNALERFC

COMMON/PROP/KSL. Q. ALPHSL. L. PI. ALPHS
COMMON/EQ/X

ALPHSR = ALPHSL*"‘O.5D0
PISR = (PI*"O.SDO)/2.D0

XX 8 X*X
X2 8 ZWX
X2ALPH 8 XZ‘ALPHSR
XALPH = X‘ALPHSR
XXALPH = XX‘ALPHSL

BETA = RAD/((4.0DO*ALPI-IS"TIME)“O.5D0)

NEWBETA = RAD/((4.0D0*ALPHS‘(TIME - TIMEC))“O.5OD0)
BETAS = BETA‘BE‘T'A

BET A2 = 2.0D0‘BET A

NEWBETS = NEWBET A‘NEWBET A

NEWBETZ = ZOWNEWBET A

DER] = Q'EXP(-BETAS)/(2.0D0*BETAS)
DER2 = -RAD/(4.0D0*(ALPHS“0.50D0)"('T'IME“(3.0D0/2.0DO)))
DER3 = Q*EXP(-NEWBETS)I(2.0D0*NEWBETS)

154

DER4 = --RAD/(4.ODO"(ALPI'IS"‘“(I.50D0)"'«TIME-THMEC)""I
+ (3.0D0/2.0D0)))
DER5 = ((-2.0IX)*BETAS‘AIPHSL‘EW-AIPHSL‘BETAS) -
+ EXP(-ALPHSL*BETAS))/(2.0DO"BETAS*ALPHSR)
+ + (ALPHSR‘EXH-ALPHSL"BETAS)»/(EXP(-ALPHSL
+ ‘XX)/(2.0D0*X*ALPHSR) - (PISR’ERFC(ALPHSR"X)))
DER6 = ((-2.0D0*NEWBETS"AIPHSL‘E)(P(-ALPHSL*NEWBETS) -
+ EXH-ALPHSL‘NEWBETSD/(Z.ODO*NEWBETS*ALPHSR)
+ + (ALPHSR‘EXPG-ALPHSL"NEWBETS)))/(EXP(-ALPHSL
+ *XX)/(2.0DO’X"ALPHSR) - (PISR‘ERFC(ALPHSR"X)))
C
IF (BETA .LT. X) THEN
FUNCLZA = DER1*DER2
ELSE
IF (BETA .GT. X) THEN
FUNCLZA = DER5'DER2
ENDIF
ENDIF
IF (NEWBETA .LT. X) THEN
FUNCLZB = DER3’DER4
ELSE
IF (NEWBET A .GT. X) THEN
FUNCLZB = DER6‘DER4
ENDIF
ENDIF
FUNCL2 = FUNCLZA - FUNCLZB
RETURN
END

SUBROUTTNE MODELZCTTMEC. TIMEC2)

DOUBLE PRECISION X. RAD. RAD2. TIMEC. TIMEC2. TIME2.
+ TIMEZB, ZBRENTZ. ZBRENTB. ERFC. TMIN
COMMON/EQ/X
COMMON/RD/RAD. RAD2
COMMON/COUNT/M
COMMON/TEMP/T'MIN

EXTERNAL ERFC. ZBRENT‘Z. ZBRENT3

CALL MODEL3(TIMEC. TMIN)

MAX = 50
DO M = 1. MMAX
IF (M .EQ. 1) THEN
TIMECZ = ZBRENT3(1.0D-4. 2.0D0. 1.0D-6. TIME2. RAD2)
ELSE
IF (M .GE. 2) THEN
TIMEZB=ZBRENT2(TIMEC2+1.0D-9,TTMEC2+0.02DO.1.0D-6.RAD2.TTMECZ)

155

TIMEC2 = ZBRENT3(1.0D4. time2B - 1.0D-6. 1.0D-6.TIME2B. RAD2)
ENDIF
ENDIF
IF (ABS(TIMEZB . TIME2) .LE. TDD-3) THEN
RETURN
ELSE
TIME2 = TIMEZB
ENDIF
IF (M .EQ. MMAX) THEN
WRITE(‘.")’NUMBER OF ITERATIONS OF M EXCEEDED’
ENDIF
ENDDO
RETURN
END

DOUBLE PRECISION FUNCTION ZBRENT3(X1. X2. TOL. TIME2. RAD2)
PARAMETER(ITMAX=1W. EPS= 3.0E-8)
DOUBLE PRECISION A. B. C. D. E. FA. FB. FC
DOUBLE PRECISION TOL]. TOL. X1. X2. XM

DOUBLE PRECISION P. Q. R. S. FUNCL3. X. TIME2. RAD2. TMIN

COMMON/EQ/X
COMMON/TEMP/TMIN
COMMON/COUNT/M
EXTERNAL FUNCL3

=X1
B=X2
FA=FUNCL3(A. TIME2. RAD2)
FB=FUNCL3(B. TIME2. RAD2)

IF(PB‘FA .GT. 0.000) PAUSE ’ROOT MUST BE BRACKETED FOR
+ ZBRENTS.’
FC=FB
DO 15 I'TER=1.1'1'MAX
IF(FB‘FC .GT. 0.0D0) THEN
C=A
FC=FA
D=B-A
E=D
ENDIF
[E(ABS(FC) .LT. ABS(FB)) THEN
A=B
B=C
C-A
FA=FB
FB=FC
FC=FA
ENDIF

156

TOLl=2.0DO*EPS‘ABS(B)-t-O.5H)"TOL
XM=0.5D0*(C-B)
1F(ABS(XM) .LE. TOL] .OR. FB .EQ. 0.0D0)THEN
ZBRENT3=B
RETURN
ENDIF

IF(ABS(E) .GE. TOL] .AND. ABS(FA) .GT. ABS(FB)) THEN
S=FBIFA
IF(A .EQ. C)T'I'IEN
P=2.0D0"'XM"S
Q=1.0D0 - S
ELSE
Q=FAIFC
R=FB/FC
P=S‘(2.0D0"XM*Q*(Q-R) - (B-A)"(R-1.0D0))
Q=(Q-1 .ODO)*(R-1 .ODO)*(S- I .0113)
ENDIF
IF(P .GT. 0) Q = -Q
P=ABS(P)
IF(ZDWP .LT. MIIVI(3.0WXM*Q-ABS(TOL1‘Q),ABS(E*Q)))THEN
E=D
D=P/Q
ELSE
D=XM
E=D
ENDIF
ELSE
D=XM
E=D
ENDIF
A=B
FA=FB
E(ABS(D) .GT. TOLDTI'IEN
B=B+D
ELSE
B=B+SIGNCTOLI.XM)
ENDIF
FB=FUNCL3(B. TIME2. RAD2)
CONTINUE

PAUSE ‘ZBRENT 3 EXCFEDING MAXIMUM ITERATIONS.’
ZBRENT3=B
RETURN
END
DOUBLE PRECISION FUNCTTON FUNCL3(TIMECZ. TIME2. RAD2)

DOUBLE PRECISION KSL. Q. ALPHSL. L. X. PI. ERFC. ALPHS. TIME2.

157

+ BETA. BETAS. BETA2. NEWBETA. NEWBETS. NEWBETZ. ETA]. ETAZ.
+ XX. X2. X2ALPH. XALPH. XXALPH. ALPHSR. PISR. RAD2. TMIN.
+ TIMECZ

COMMONEQ/X

COMMON/TEMP/TMIN

COMMON/PROP/KSL. Q. ALPHSL. L. P]. ALPHS
COMMON/COUNT/M

ALPHSR = ALPHSL“0.5D0

PISR I: (PI"0.5D0)/2.D0

XX = X‘X

X2 = 2.0DO*X

X2ALPH = XZ‘ALPHSR

XALPH = X‘ALPHSR
XXALPH = XX‘ALPHSL

IF (M .EQ. 1) THEN
TIME2 = 1.050D0*TIMEC2
ENDIF
BETA = RAD2/((4.0D0‘AIPHS‘TIME2)“0.5D0)
NEWBETA a RAD2/((4.0D0‘ALPHS‘TTTME2 - TIMECZ))"0.50D0)
BETAS = BETA*BETA
BET A2 = 2.0DO*BETA
NEWBETS = NEWBET A‘NEWBET A
NEWBETZ = 2.0DO"NEWBETA

IF (BET A.LT.X)THEN
ETA] = 1.0D0 - Q"(EXP(-BETAS)/BET A2 - EXP(-XX)/X2 -
+ PISR*(ERFC(BET A) - ERFC(X)»
ELSE
IF(BETA.GT.X)TT‘IEN
ETA] = (EXH-ALPHSL‘BETASMBETAZ‘ALPHSR) -PISR*
+ ERFC(AIPHSR*BETA))/(EXP(—XXALPH)IX2ALPH
+ -PISR"'ERFC(XALPH))
ELSE
ETA] a 1.0D0
ENDIF
ENDIF
IF (NEWBETA .LT. X) THEN
ETAZ = 1.0D0 - Q‘(EXP(-NEWBETS)/NEWBET2 - EXP(-XX)IX2 -
+ PISR*(ERFC(NEWBET A) - ERFC(X)»
ELSE
IF (NEWBETA .GT. X) THEN
ETAZ :3 (EXP(-ALPHSL“‘NEWBETS)I(NEWBET2*ALPHSR)
+ -PISR*ERFC(AIPHSR*NEWBETA))I
+ (EXP(-XXALPH)/X2AIJ’H -PISR“ERFC(XAIJ’H))
ELSE
ETA2 = 1.0D0

C

DO

0000

158

ENDIF

ENDIF

FUNCL3 a: TMIN- (ETAI - ETA2)
RETURN
END

SUBROUTINE MODELMTIMEC. TMIN)

IMPLICIT REAL'S (A-H.O-Z)

DIMENSION 'r(3500.5).v(3500).srcz(3500).B(5).2(5).
+A(5).BS(5).VINV(5.5).EXTRA(20)

DIMENSION P(5.5).PS(5.5)

DOUBLE PRECISION KSL. Q. L. PI, ALPHSL. ALPHS
DOUBLE PRECISION ERFC. ZBRENT]
DOUBLE PRECISION ALPHSR. PISR. X. XX. X2. X2ALPH.
XALPH. XXALPH. BET A. BET AS. BETA2. ETA. NEWBETA.
NEWBETS. NEWBETZ, TIME. TIMEC. RAD. ETA]. ETA2.
RAD2

COMMON SIGZ.T.2.BS.I.ETA.PS.P.B.A.Y.MODL.VINV.NP
+.EXTRA
(DMMONIMOD/AAJL

COMMON/PROP/KSL. Q. ALPHSL. L. PI. ALPHS
EXTERNAL ERFC. ZBRENT]

KSL 8 LOW
ALPHSL = 1.0D0
L = -1(X).0D0

Q = -1.0D0

ALPHS = 1.0D0
RAD = T(I,2)
RAD2 = T(I.1)

PI = DACOS(-1.D0)
TTMEC = BS(1)

ALPHSR = ALPHSL**O.5D0

PISR = (Pl*‘0.5D0)/2.D0
X IS THE CALCULATED VALUE FOR LAMBDA
ZBRENTisaNfindhgsubmufineusedmsolvethehanscendanalequafionforme
freezing from location. See Numerical Recipes by Press et al.. Cambridge University
Press. New York. New York. 1986.
X a ZBRENTl(l.0D-4.2.D0.0.001D0)

XX = X‘X

X2 = 2.DO*X

X2ALPH = X2*ALPHSR

XALPH = X‘ALPl-ISR

XXALPH = XX*ALPi-ISL

TIME = ZBRENT2CT'IMEC+1.0D-9.TIMEC + 0.010D0.1.0D-6. RAD. TTMEC)

C

000

159

BETA = RAD/((4.0DO*ALPHS"'T'IME)“O.SIXI)
NEWBETA = RAD/((4.0D0‘ALPHS‘1TIME - TIMEC»“O.50D0)

BETAS = BETA’BETA

BET A2 = 2.D0"BETA

NEWBETS = NEWBETA‘NEWBETA

NEWBETZ = 2.0D0‘NEWBET A
IF (BET A.LT.X)THEN

ETA] = 1.0D0 - Q*(EXP(-BET AS)/BET A2 - EXP(-XX)/X2 -
ELSE PISR‘(ERFC(BETA) - ERFC(X)»

IF(BETA.GT.X)THEN
ETA] = (EXP(-ALPHSL*BETAS)I(BETA2*ALPHSR) -PISR*
ERFC(AIPHSR‘BETADKEXH-XXALPHMXZALPH

-PTSR"'ERFC(XALPH))
ELSE
ETA] = 1.0D0

ENDIF
ENDIF
IF (1 .LE. EXTRA(1)) THEN
ET A2 = 0.0D0
ELSE

IF (NEWBETA .LT. X) THEN
ETA3 = EXH-NEWBETS)/NEWBET2
ETA4 = EXP(-XX)/X2
ETA5 a ERFC(NEWBETA)
ETA6 = ERFC(X)
ETA7 = 1 - Q‘(ETA3 - ETA4 - PISR“(ETA5 - ETA6))
ETA2 = ETA7
ELSE
IF (NEWBET A .GT. X) THEN
ETA2 = (EXH-ALPHSL‘NEWBETSWNEWBETZ‘AIPHSR)
-PISR*ERFC(ALPHSR*NEWBETA))I
(EXP(-XXALPH)/X2ALPH -PISR*ERFC(XALPH))
ELSE
ETA2 a 1.000
ENDIF
ENDIF
ENDIF
ETA = ETAl-ETAZ
TMIN = ETA
RETURN
END

160
ThisfilerepresentstheinputfileforusewithNLlNAfORinthedeterminationofrcwim

prior information from two difierent radius locations. The first two entries in the second column
represemthetreannent times fromtheprevious procedures,whilethefirsttwo entriesinthefourth
cohmnmpmsemmetwodifiemmmdiuslmafionwithmndomenom.1hisfifmwumnisdn
original radius values. with random em.

12 1 2 200 1 1
.125000005+00

\OGQO‘MAwN—I

10
11
12
]
2

. rosooooora+oo
225000me
.23941400E+01
.239414(X)E+Ol
.23941400E+01
.23941400E+01
.2 394 “(DE-+01
.99467(X)0E+(X)
.99467000E-t-(X)
.994670(DE+CX)
.994670(DE+(X)
.99467(XX)E+(X)

AWE-03
.100000005-03
. rooooomra-oz
. IWE-OQ
.rooooooosoz
AWE-02
. ](XDOOME—OZ
. “HOWE-02
.1WE-02
.rooooooor-zm
AWE-02
.lmme-OZ

.10167473E+(X)
.10745657E-t-(X)
.93793248E-01
.1m74897E-t-(X)
.10728538E+(X)
.rmssmmoo
.82217395E-01
.1W87963F.+00
.86525441E-01
.95035257E-01
.91867172E-01
.92593108E-01

.lm98 108m
.971 19240E-01
.95522674E-01
.99939520E-0]
.10820452E+(X)
.99252245E-Ol
.84363144E-01
.168289'79E+m
. 15766006E+m
.16359131E+(X)
.14363310E+00
.16108979E+m

BIBLIOGRAPHY

Augustynowicz. SD. and Gage. AA. 1985. 'Temperatme and Cooling Rate Variations During
Cryosurgical Probe Testing." International Journal of Refrigeration. Vol. 8. pp. 198-208.

Beck. JV. 1991. NLINA.FOR. Fortran Computer Program. Department of Mechanical
Engineering. Michigan State University. ’

Beck. J.V.. and Arnold. 10., 1977. Parameter Estimation in Engineering and Science. John
Wiley & Sons. New York.

Budman. B.. Dayan, J.. and Shitzer. A.. 1990. "Controlled Freezing and Thawing of a Non-Ideal
Solution with Application to Cryosurgery.” Single and Multiphase Convective Heat
Transfer. ASME. Vol. 146. pp. 55-59.

Comini. G.. and Del Guidice. S.. 1976. ”Thermal Aspects of Cryosurgery.” Journal of Heat
Transfer. Vol. 98. no. 4. pp. 543-549.

Cooper, TE. and Petrovic. W.K.. 1974. ”An Experimental Investigation of the Temperature Field
Produced by a Cryosurgical Cannula.” Joumal of Heat Transfer. VoL 96. no. 3. pp. 415-
420.

Cooper. TE. and Trezek. 6.1.. 1970. ”Analytical Prediction of the Temperature Field Emanating
from a Cryogenic Surgical Cannula.” Cryobiology. Vol. 7. no. 2-3. pp. 79-93.

Cooper. TE. and Trezek. 6.1.. 1972. "On the Freezing of Tissue.” Journal of Heat Transfer. Vol.
94. no. 2. pp. 251-253.

Filippov. Y.P.. and Vasil'kov. A.P.. 1978. "Tissue Temperature Simulation in Cryosurgery."
Journal of Engineering Physics. Vol. 36. no. 6. pp. 725-729.

Gage. AA. 1992. "Cryosurgery in the Treatment of Cancer.” SURGERY. Gynecology &
Obstetrics. Vol. 174. pp. 73-92.

Gilbert. J.C.. Wk. G.M.. Haddick. WK. and Rubin8ky. B.. 1984. "The Use of Ultrasound

ImagingmeommnngCryosmgery."Pmceedingsofme6mAmualIEEE/Engineenng
in Medicine and Biology Society Conference. New York. pp. 107-111.

Gilben. I.C.. Rubinsky. B.. and Onik. GM. 1985. ”Solid-Liquid Interface Monitoring with

Ultrasound During Cryosurgery.” Submitted to the American Society of Mechanical
Engineers Winter Annual Meeting. paper no. 85.

161

162

Hayes. L.J.. and Diller. KR. 1982. ”A Finite Element Model for Phase Grange Heat Transfer in
11 Composite Tissue with Blood Perfusion.” Modeling and Simulation. Proceedings of the
13th Annual Pittsburgh Conference. Vol. 13. pp. 59-63.

Hrycak. R. Levy. MJ.. and Wilchins. SJ.. 1975. ”Cryosurgery of Lesions Through Contact
Freezing and Estimates of Penetration Times.” Submitted to the American Society of
Mechanical Engineers Winter Annual Meeting. paper no. 75.

Keanini. R.G., and Rubinsky. B.. 1992. "Optimization of Multiprobe G'yosurgery.” Journal of
Heat Transfer. Vol. 114. no. 4. pp. 796-801.

Ozisik. M.N.. 1980. Heat Conduction. Jolm Wiley & Sam. New York.

Paterson. S.. 1952. "Propagation of a Boundary of Fusion.” Proceedings of the Glasgow
Mathematical Association. Vol. 1. pp. 42-47.

Press. W.H.. Vetterling. W.T.. Teukolsky. and S.A.. Flannery. B.P.. 1986. Numerical Recipes,
Cambridge University Press. New York.

Rubinsky. B.. 1986a. "Cryosurgery Imaging with Ultrasound.” Mechanical Engineering. Vol. 108.
pp. 48-52.

Rubinsky. B.. 1986b. ”Recent Advances in Cryopreservation of Biological Organs and in
Cryosurgery." Heat Transfer. Proceedings of the 8th International Heat Transfer
Conference. Vol. 1. pp. 307-316.

Rubinsky. B.. and Eto. K.. 1989. ”Heat Transfer with Phase Transition in Biological Materials."
Cryo-Letters. Vol. 10. no. 3. pp. 153-168.

Rubinsky. B.. and Shitzer. A.. 1976. ”Analysis of a Stefan-Like Problem in a Biological Tissue
Around a Cryosurgical Probe." Journal of Heat Transfer. paper no. 76. pp. 514—519.

Savic. M.. 1984. ”Microprocessor Enhanced Accuracy of Cryodestruction.” Proceedings of the 6th
Annual [BEE/Engineering in Medicine and Biology Society Conference. New York. pp.
79-81.

Scott. E.P.. 1993. "An Analytical Solution and Sensitivity Study of Sublimation-Dehydration
Within a Porous Medium with Vohunetric Heating.” Submitted for publication to the
Journal of Heat Trrmsfer. January. 1993.

Wanen. R.P.. Bingham. PE. and Carpenter. JD. 1974. "Heat Flow in Living Tissue During
Cryosurgery." Submitted to the American Society of Mechanical Engineers Winter Annual
Meeting. paper no. 74.

"Iiiiitiii’ii7tii7fliis