OF CATERPILLARS AND MODELS by Richard T e r r e l l Steinhoff B.Sc, Colorado State University, 1970 THESIS SUBMITTED IH PARTIAL FULFILHENT OF THE REQUIREHENTS FOB THE DEGREE OF MASTER OF SCIENCE i n the Department of COMPUTER SCIENCE We accept t h i s thesis as conforming to the required standard The University of B r i t i s h Columbia October 1975 R I G H T S OF P U B L I C A T I O N AND L 0 AN In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I - agree that the Library shall make i t freely available for reference and study.. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his representative.. It i s understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Computer Science The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date October 9, 1975 i i i A B S T R A C T An approach to modeling in general, i t s background, and i t s implementation can be described by examining i t s application to a particular system. . The Western Tent Caterpillar i s a productive subject of study because there i s extensive data available for determining the relationships within the system. . Individual differences between caterpillars and the weather are the dominant factors involved. They are incorporated into a set of ordinary differential equations that express the rate of change in the caterpillar population over time., These equations are numerically integrated to calculate the population at some terminal time. The results over a six year period agree reasonably well with . observed f i e l d data. . The three dimensional plot of the population gives a visual appreciation for the distribution over time. , i v T A B L E O F C O N T E N T S 1. , INTRODUCTION .......................................... 1 2. THE BACKGROUND - A Blending Of The Sciences ....... , . . .4 2.1 The Levels Of Application ........................4 2.2 Operations Research ..............................5 2 .3 Biology ..........................................7 2.4 Mathematics And Computer Science .................9 2 i5 The History Of The Tent C a t e r p i l l a r Model ....... 15 3. THE SET - The Tools For This P a r t i c u l a r Problem . . . . . . 1 8 3.1 S p a t i a l Boundaries .............................. 19 3.2 Components Of The System ........................ 19 3.3 State And Driving Variables 20 3.4 The Independent Variable ........................ 21 3.5 Continuous Or Discrete .......................... 21 3.6 Stochastic Elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 4 4. THE SCRIPT - In Retrospect . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 6 4.1 Determination Of The Major State Variables . . . . . . 2 7 4.2 Description Of The Model . . ...... . . . . . . . . . . . . . . . . 3 2 4 .3 Determination Of The Parameter Values . . . . . . . . . . . 3 4 5., VALIDATION ...........................................40 6. , CONCLUSIONS . . . . . . . . . . . . 4 9 7. , BIBLIOGRAPHY .........................................51 8. . APPENDICES ........................................... 52 8.1 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2 8.2 Newton-Raphson Routines ......................... 53 8.3 Constrained Targeting Routines . . . . . . . . . . . . . . . . . . 5 9 8.4 C a t e r p i l l a r Model Code ............ ,125 V L I S T OF F I G U R E S I - State Variables In Average Heather ...................38 II - Derivatives Of State Variables In Average Heather .,,39 III - Comparison Of Real And Simulated Populations .......1*2 IV - Heather And The Simulated Caterpillar Population .,.,44 V - Caterpillar Distribution In Mild Heather .............45 VI - Caterpillar Distribution In Harsh Heather ...,,..,-,,46 VII - Caterpillar Distribution In Mild Heather 47 VIII - Caterpillar Distribution In Harsh Heather .,,.,,.,,48 v i A C K N O W L E D G E M E N T S Anyone who has written a masters th e s i s knows the large debt I owe to to my . advisors and friends for guiding and encouraging me through t h i s endeavor. Roger Hart and Ilan Vertinsky . have been helpful beyond anything I had expected. John Coulthard and Mark Dionne supported my graphic display e f f o r t s , j, Without them and the f a c i l i t i e s a vailable at the University of B r i t i s h Columbia Computing Center, t h i s would not have been possible. B i l l Hebb's documentation program has made the laborious part of writing t h i s just bearable. Since Ginny and I are one i n a l l we do as husband and wife, she i s the unsigned co-author of t h i s thesis and f u l l y deserves her half of the c r e d i t ! The data used i n building the model was due to the f i e l d work of B i l l Wellington and organization of B i l l Thompson., Together with Peter Cameron and Abe Lansberg, they kept my biology s t r a i g h t , my nose i n the wind, and my eyes on the horizon. / 1 Is. INTRODUCTION Man i s constantly trying to formalize the rather nebulous and complex world around him so that he can understand i t , and manipulate i t . This formalization, i n the broadest sense of the word, may take many forms: everything from poetry to detailed mathematical proofs. One of the intermediate forms i s modeling - the a r t of defining a transformation from one space in t o another one, such that an en t i t y i n the o r i g i n a l space has, i n some sense, an "equivalent" entity i n the second space which lends i t s e l f more e a s i l y to rigorous examination. By "equivalent", I mean "responds to a given set of s t i m u l i i n a s i m i l a r manner - both quantitatively and q u a l i t a t i v e l y . " Mathematical modeling means that the transformation i s from the physical world into a mathematical representation of i t . . This requires an e x p l i c i t statement of the rela t i o n s h i p s that e x i s t within the e n t i t y being modeled. My p a r t i c u l a r application of modeling i s the western tent c a t e r p i l l a r . This was f o r a variety of reasons. One of the deciding reasons was that that area of application had to be meaningful to the world around me. Research In some areas would not only be useless and redundant, but t o t a l l y incapable of being validated. , When I found that there was some in t e r e s t , and indeed some exc i t i n g progress i n modeling tent c a t e r p i l l a r s , that became the area of study. There were several other advantages i n making t h i s choice. In b i o l o g i c a l modeling one of the biggest hindrances i s the lack of data on some of the important relationships i n the system - missing 2 jigsaw pieces, i f jou w i l l . This lack i s not due to any b i o l o g i c a l oversight, but rather to the f a c t that such data, before modeling came of age, was too massive and complex to be assimilated by an i n d i v i d u a l where as i t i s e a s i l y "digested" by including i t i n the description of the o v e r a l l model., The advantage with the tent c a t e r p i l l a r s was that H.G. Wellington [ 6] has devoted a great deal of time i n establishing the necessary r e l a t i o n s h i p s needed to describe the tent c a t e r p i l l a r s . , This thesis w i l l have three major thrusts., One feature of the thesis as a whole i s that i t w i l l present a step by step approach to b i o l o g i c a l modeling from the computer s c i e n t i s t s viewpoint. , I w i l l indicate i n the f i r s t chapters the general d i s c i p l i n e s and tools needed to attack t h i s area of research with . s p e c i a l attention given to d i g i t a l simulation. , I was fortunate to be i n a department with some f l e x i b i l i t y and a v i s i o n for the wide range of applications of modeling. , By . e n r o l l i n g i n courses i n other f i e l d s , I have come to some conclusions as to what others would f i n d h e l p f u l i f they were to study bio-simulation. I therefore submit these i n the hope that others may learn by t h i s experience., although - 1 w i l l be discussing p a r t i c u l a r s i n the l a t t e r chapters, I hope that the approaches outlined there w i l l be of some use i n other applications., The second goal i s to r e l a t e the various types of modeling and itemize t h e i r strengths and weaknesses to a variety of problems. , From that I w i l l present the important 3 c h a r a c t e r i s t i c s of tent c a t e r p i l l a r s and describe the type of modeling I chose for them.. F i n a l l y , I w i l l discuss my model in d e t a i l , with s p e c i a l attention given to how the variables and constants of the system were derived and a presentation of the evidence supporting .the hypothesis that t h i s i s a r e a l i s t i c model. In retrospect, i t ' s i n t e r e s t i n g that c a t e r p i l l a r s should be the p a r t i c u l a r application of modeling that was most promising a f t e r examining a number of possible areas. The metamorphosis that the c a t e r p i l l a r undergoes i n changing into a b u t t e r f l y or moth has been the analogy used f o r r e b i r t h . Jesus t o l d a f i r s t century scholar, "you must be born again" (JOHN 3:3) and to the scholar that was completely incomprehensible, just as the c a t e r p i l l a r cannot f l y u n t i l i t becomes a completely new . creature and transcends the l i m i t a t i o n s i t had as a c a t e r p i l l a r . Wellington expressed that need i n research "at the border of the unknown, one must consciously s t r i v e t o escape from the mesh of former frames of reference, and to remain outside the generally accepted range of opinion concerning one's problem..." The c a t e r p i l l a r has the key to a s p i r i t u a l truth which . i s completely inapproachable by s c i e n t i f i c i n v estigation, , Unless we escape our t r a d i t i o n a l " s c i e n t i f i c " frame.of reference when we come to explore the unknown side of our being - our s p i r i t s - we w i l l never possess our wings to carry us above our reason! 4 Is. , THE BACKGROUND - A BLENDING OF THE SCIEHCES Mathematical modeling may be approached on one of a variety of levels and incorporate a vide spectrum of disciplines., In this chapter I will discuss the levels and the variety . of s k i l l s from both the biological and physical sciences necessary to deal with the problems presented at the various levels. , The history of the tent caterpillar study i s also included., 2.1 THE LEVELS OF APPLICATIOH In the attempt to build an overview of the structure these levels define, l e t me indicate the two extremes.. The most detailed level i s the study of the relationships between component parts of the system. In a biological environment (rather than, for example, a bio-chemical - one), this would usually mean the guantitative effect of one biological element upon another., For example, the mountain lion population in a given area effects the deer population that area. This effect could be expressed as either a death rate or as a discrete event, that being the death of a deer., On this detailed level, the study would center on a particular aspect of the effect the mountain lion has on a deer. For example, i t might focus on the effects on i t s feeding habits and model that relationship. This i s in contrast with modeling the deer and mountain lion population as a whole. On the other extreme, the macro level, the model as a whole could be used to make management decisions about the system, based upon certain constraints and goals to be 5 s a t i s f i e d . , A f i s h and game department might decide upon the number of game permits to be issued by finding an "optimal*! strategy to harvest the deer, y This strategy , could be formulated and tested i n a model, , The intermediate l e v e l s usually t r y . to study o v e r a l l trends i n the system by reaching a f a i r l y tenable description of the system and then a l t e r i n g the environment i n which the system i s operating. , A p a r t i c u l a r model f o r deer may be subjected to a variety . of catastrophic conditions (such as fo r e s t f i r e s , epidemics, etc) and the recovery process examined, , The following chapter w i l l discuss t h i s intermediate l e v e l rand the d e t a i l e d l e v e l , since they were the areas of concern f o r t h i s t h e s i s . . Let me b r i e f l y describe the l e a s t detailed l e v e l and the kinds of research i n that area.. 2.2 OPERATIONS RESEARCH This macro l e v e l i s known most frequently as operations research. , I t was born during World War I I as a r e s u l t of the problem of where to place a f i n i t e number of radar antennae such that the enemy airplanes could be most e f f i c i e n t l y detected.,• Operations research i s usually employed f o r the expressed purpose of making management type decisions.,- As a ru l e i t does not require an in-depth understanding of the model i t s e l f , but depends upon i t to react i n the same way that the r e a l l i f e system would to external s t i m u l i . r From these reactions i t forms a set of p o l i c i e s to be followed that w i l l cause, the system to behave i n a specified.manner. These 6 p o l i c i e s might e i t h e r be of the yes-or-no variety, or they might be of the quantitative type - harvest so many in d i v i d u a l s . , A further complicating f a c t o r i s the discrete nature of some decisions; one cannot harvest 1.5 i n d i v i d u a l s . The p o l i c i e s may take the from of simple constant terms (harvest 400 deer) or may depend upon the current state of the system (how many deer are present, how many mountain l i o n s e t c ) . These dependent p o l i c i e s are usually expressed as a weighted sum of the dependent variables or as a table of points which represent the given policy plotted against the possible state of the system. Anyone interested i n t h i s aspect of modeling should have some background i n operations research and f a m i l i a r i t y with l i n e a r and non-linear programming and optimization theory. The material i n t h i s area describes a tones that#which# r e a l world) to optimize some performance function subject to a set of constraints, e i t h e r on variables i n the system or on the control parameters. This usually requires the construction of a Jacobian (a matrix of p a r t i a l derivatives of the constraints and performance index r e l a t i v e to these control parameters) and a set of operations on i t . Hence, a strong background i n 1 Courses at the University of B r i t i s h Columbia that would be useful for those interested i n Operations research are; COMMERCE 310,410,411,581,582,584,586,590>681,682; COHP0TER SCIEHCE 406; ELECTRICAL ENGINEERING 467,485,555,557,568,570; FORESTRY 331,559; MATHEMATICS 140,151,156,221,321,362. 7 matrix theory and vector analysis i s also d e s i r a b l e 1 . 2.3 BIOLOGY There are some general background areas i n biology that are h e l p f u l i n modeling ecosystems. „ Obviously, the p a r t i c u l a r biology f o r a given system must be f u l l y understood before the important components of the system can be i s o l a t e d and t h e i r i n t e r a c t i o n can be formalized. y I t i s v i t a l that the computer s c i e n t i s t confers v i t h a b i o l o g i s t who i s knowledgeable i n the problem. Indeed, one of the side e f f e c t s of modeling that w i l l be. discussed l a t t e r i s that new d i r e c t i o n s of study are needed to complete the unknown relationships discovered during the formalization of the model.,, Formalization often points out unsubstantiated assumptions i n any d i s c i p l i n e . The general b i o l o g i c a l background that would be useful to a computer s c i e n t i s t i n the understanding an ecosystem, would s t a r t with a basic understanding of botany and zoology, and then continue with advanced work which focuses on s p e c i f i c ecosystems., Special attention should be given to the transfers of energy and mass between components i n the system and chemical cycles. This exposure to the complicated variety of relationships - along with . a f e e l i n g f o r the r e l a t i v e importance of them i s the r e a l basis f o r building a b i o l o g i c a l 8 model** At this point - the formalization of the "biology" in the model - one has to confront the rather obvious question of the ab i l i t y of a computer scientist to deal with some issues that are biological in nature. / Science i s a spectrum of knowledge in which each of the basic disciplines overlaps. ,• In the purest form, a mathematician or computer scientist i s a craftsman - turning out tools for understanding and manipulating the world around us., As soon as he becomes involved with a particular application, he begins to become a biologist, engineer, or social scientist.„ It has been argued that there i s no "practical" need for the purist, but that neglects the contribution of the tool i t s e l f , i t s study, improvement and upkeep. , So the computer scientist and mathematician can give to a biologist some well honed tools numerical analysis, mathematical model construction, and simulation languages which incorporate the whole gamut of theoretical design and programming technigue. It must be stressed that this i s a far cry from the s k i l l s traditionally associated with a computer programmer, in the same way that a mathematician has grown past the multiplication table stage of early c i v i l i z a t i o n . The biologist and computer scientist. * Courses at the University of British Columbia that would be informative for computer scientists who want to apply their modeling to biology would be: APPLIED SCIENCE 250; BIOLOGY 301,321,322; .FORESTRY 405,409,505; ZOOLOGY 401,421,527. 9 then, need a continual dialogue as they approach the problem of bio-simulation. ., There i s one additional advantage by having a computer scientist approach the biological problem., That i s : he does not have any preconceived notions or conceptual frameworks to bias his view of the biological "facts". , For example. Hay [3] has written one of the most comprehensive and accepted summaries of the state of the art biological modeling, and his background i s in ele c t r i c a l engineering! The above comments should not be misconstrued to mean that the biological understanding i s not a requirement for successful modeling, but merely points out the need for cooperation between computer scientists and biologists in the f i e l d of bio-simulation. , 2.4 MATHEMATICS AND COMPUTER SCIENCE Something should be said about the contributions to modeling and computer simulation by mathematics and computer science., I have divided the areas of effort into numerical analysis, general algorithms, theory of simulation, and language design. The s k i l l s mentioned in these areas should be regarded as minimal recommendations for any computer scientist who aspires to bio-simulation 1. 1 Specific courses at the University of British Columbia that should be the basis for a curriculum in bio-simulation are: COMPUTER SCIENCE 115,215,315,350,405,505,535; ELECTRICAL ENGINEERING 591; MATHEMATICS 155,255,256,301. 10 Numerical analysis has been one of the most important aids to successful modeling. For the continuous system, numerical integration and d i f f e r e n t i a t i o n techniques have been improved, so that they are more accurate and e f f i c i e n t . They also have the a b i l i t y to dynamically change the integration step s i z e depending upon the non-linearity of the rate equations. The most advanced techniques i n t h i s area are Predictor-Corrector and Runge-Kutta methods with eigenvalue step s i z e modification. Advanced work has been done i n handling " s t i f f " systems of d i f f e n t i a l equations. , The question of whether a system i s " s t i f f " or not can generally be l e f t u n t i l one actually begins to simulate the system. ., I f i t i s a s t i f f system, there should not be any doubt, once the integration has begun, because the state variables w i l l become very large i n absolute value. Since eigenvalues are used to a l t e r the step s i z e s during numerical integration, a set of eigenvalue routines f o r general and s p e c i a l matrices have been developed. Often i t i s necessary . to incorporate table orientated data i n t o continuous system simulation programs. Numerical analysts have developed curve f i t t i n g to approximate the values i n the table with l e a s t squares algorithms and spl i n e s . , l i n e a r equations that a r i s e i n simulation problems require methods to solve that s p e c i a l type of system i f i t can not be solved a n a l y t i c a l l y , so mathematicians developed proqramming and integer programming* State of the art i n t h i s area i s the Revised Simplex Algorithm by Danzig. The c l a s s of general algorithms or program techniques 11 that are particularly useful in simulation includes sorting, hashing, and data handling. Hashing i s the process of transferring large set of entities into a small set of classes such that the number of classes i s in some sense optimal, and that roughly the same number of entities are l i k e l y to f a l l i n each class. Once a hashing algorithm has been decided upon for a given set of entities, one can search the set for a particular entity more efficiently by f i r s t going to the appropriate class. Hashing vould be useful in discrete simulation problems that reguire one to find a given individual amount a large number of similar individuals. For example, in modeling a employment agency, one might want to find a particular job applicant in a f i l e of several thousand applicants. Sorting i s the process of making an ordered l i s t of entities., This could be done by a variety of methods such as binary search, bubble sort, and sort/merge. Sorting i s useful in block diagram simulations where the blocks must be sorted before they can used to direct the simulation. .. Data representation i s one of the areas of most interest in computer science., Several special types have already been examined and used in this thesis. Trees, l i s t s , and multiple-linked structures are a l l very useful in simulation. The linked l i s t used to connect occupied grid squares in the caterpillar model saved unnecessary integration time. In the purest form, simulation has been studied from several theoretical aspects. Stability i s a crucial question in modeling.„ The mathematician i s able to determine the 12 s t a b i l i t y of some systems a n a l y t i c a l l y and the computer s c i e n t i s t has developed steady state techniques to i s o l a t e transients i n the system and determine t h e i r ultimate e f f e c t . Deadlock i s another area of study., This i s the state where two or more processes freeze uncompleted because neither can proceed u n t i l the other process i s completed. As an example of t h i s , consider the t r a f f i c s i t u a t i o n where two cars want to turn l e f t just a f t e r they pass each other, one coming from the east and one coming, from the west., A l i n e of cars begins to form behind both cars., Both cars are prevented from turning by the other car's lineup, and that lineup w i l l not move u n t i l they turn* The question of deadlock i s more prevalent i n a discrete simulation, and involves a p r o b a b i l i t y rather than a absolute state of deadlock or non-deadlock.. Validation, which i s also a t h e o r e t i c a l area of computer science, w i l l be discussed i n chapter V. . The l a s t area of computer science I want to describe i s language design and implementation. That i s an important subject given the following two f a c t s : the computer i s merely a very f a s t and accurate information processing machine with a memory for storing the input data, the r e s u l t s , and the instructions for which operation to do next, a l l of which ultimately are represented as numbers; humans recognize and think i n symbols and nested processes relevant to a given problem or a c l a s s of problems. I t i s therefore important to have a language i n which the problem (which i n t h i s case i s simulation) can be e a s i l y formulated and f o r which there 13 exists an unambiguous mapping between that language and the computer*s i n t e r n a l language (numbers). A detailed history of computer languages i s out of the scope of t h i s t h e s i s . The interested reader i s referred to Jean Sammet's book [3] on programming languages and t h e i r history* In computer science, languages evolve i n a s u r v i v a l of the f i t t e s t environment. Old and unused languages are usually replaced by more powerful ones.„ There i s a ce r t a i n amount of i n e r t i a i n the system, but eventually, the poorer languages tend to die out. What follows i s a subjective judgement on the most useful . languages for simulation. No doubt i n the near future people w i l l read these l i n e s and smile at the primitive languages suggested. There are four relevant types of languages f o r simulation: general, discrete, block diagram, and continuous. The f i r s t category . are languages f o r a broad range of s c i e n t i f i c applications., ALGOL, and i t s predecessor FORTRAN, are the two most widely used. The design c r i t e r i a f o r ALGOL-68 (which i s an improved version of ALGOL) has been suggested, but the generality and complexity of the language have made i t d i f f i c u l t to implement and there are to date no complete versions available. SIHOLA-67 contains not only the a b i l i t i e s of ALGOL, but co-routines which could be used to simulate p a r a l l e l processes.. Unfortunately, i t i s not very widely used i n comparison, but as i t comes into i t s r i g h t f u l appreciation by the s c i e n t i f i c community, i t w i l l . be used more and therefore be the best choice for a general language. 14 Discrete languages would be used for discrete simulation which w i l l be contrasted with continuous simulation in chapter III* , The most: useful languages in this class are SIHSCIPT and GPSS. . GPSS i s a block orientated language, the most logical, and the easiest to learn. SIHSCBIPT i s a statement orientated language and more flexible than GPSS. / SIMULA-67 also contains discrete simulation f a c i l i t i e s . Block . orientated languages (the predecessor of the continuous system simulation languages) are used for a class of problems originally handled by an analog computer. The best of this class in common use i s DSL/90. Its predecessors were MIDAS and MIMIC. , Finally, continuous system languages in common use are DYBJAMO and CSHP, the latter being far and above the most widely used., For simple problems, these are sufficient, but there i s s t i l l a need for a more powerful, interactive continuous system simulation language. A report on the desirable features of such a language was produced in 1967 [5 ] . This CSSL report specified design objectives that included the better features of CSMP, DSL/90, GPSS, and ALGOL. It also made an attempt to extend the f l e x i b i l i t y of these languages to user defined operators and control structure. .. It promised graphic display of output and a conversational mode during the simulation. A f u l l blown implementation has not been forthcoming, but the objectives are clear and one would expect in the near future to see a third generation language: CSSL. 15 2,5 THE HISTORY OF THE TENT CATERPILLAR MODEL While i n the general overview of modeling, i t i s h e l p f u l to look . at some examples of applications of modeling. The range of modeling goes from p o l i t i c a l or s o c i a l . science and armament i n developing countries to nuclear physics and control of the f i s s i o n process i n reactors. I t s domain goes from energy cycles within the c e l l to examining the big bang theory i n space. t The c a t e r p i l l a r problem i t s e l f has some history that should be related here. The Western Tent C a t e r p i l l a r , Malacosma pluvialg (Dyar), has a population that fluctuates quite a b i t from year to year. On the Saanich Peninsula of southeastern Vancouver Island, the population varied between several thousand to several m i l l i o n i n d i v i d u a l c a t e r p i l l a r s over a period of ten years [10]., That increase and decrease over , the years became the area of study i n 1955 f o r several entomologists [6] and the subject of many subsequent papers i n population dynamics [8] [9] [ 1 ] . , The i n i t i a l question that was asked was "why don ,t the c a t e r p i l l a r s continue to increase without bound?" At the time that the question was f i r s t asked, the leading theory among those who study population dynamics was that populations were controlled by density dependent "regulators" i n the ecosystem. These regulators might be predators of the population or the food source of the population. The idea i s that as the population changed, there was a corresponding e f f e c t upon the predators or food source. I f , f o r example, predators increase by some exponential 16 function of the number of prey, the prey would eventually cause enough predators to be around so that the prey are eaten fas t e r than they can reproduce* The prey population would therefore be c o n t r o l l e d . The key word here i s dependence. Wellington [8] was just beginning to look at the tent c a t e r p i l l a r s i n 1948 with the idea of studying the e f f e c t s of weather on i t s population* The more he observed them, the more he became convinced that density dependent regulation was not the only factor i n the o s c i l l a t i o n of populations. In the f i r s t place i t i s very clear that i n d i v i d u a l animal behavior and a c t i v i t y influences i t s ultimate s u r v i v a l . , He succeeded i n showing the variation i n a c t i v i t y l e v e l within the tent c a t e r p i l l a r . population and demonstrated i t s e f f e c t on o v e r a l l colony s u r v i v a l . , Further, although predators were few and food sources were usually unlimited, the c a t e r p i l l a r population did reach a peak i n 1956 and began to decline. I t had been his conviction for some time that weather was a very large determining f a c t o r i n c a t e r p i l l a r mortality, and t h i s became more substantiated as his studies progressed. So there were two major break-throughs i n t h i s study: i n d i v i d u a l t r a i t s e f f e c t i n g o v e r a l l population growth, i n addition to weather, which i s a density independent fact o r . . Host of his observations between 1955 and 1965 exploited t h i s new knowledge and were consequently quite he l p f u l i n constructing a model which incorporated these f a c t s . , His v i s i o n for the important re l a t i o n s h i p s within the c a t e r p i l l a r population was the reason for the unusual amount of data which 17 now can be used to model the c a t e r p i l l a r ecosystem., Subsequently, P. Cameron and V. Thompson et a l [ 7 ] , began to b u i l d a model f o r the c a t e r p i l l a r population. They had the i n i t i a l task of transforming the f i e l d data on micro climate ( l o c a l climate), tree height, egg masses, and i n i t i a l populations to a form that could be rea d i l y used i n a computer program. , Their work i n that area was invaluable to t h i s thesis., Once that was done, they completed a model fo r the c a t e r p i l l a r system. , B a s i c a l l y , i t was a continuous simulation, although the rate equations were difference equations and the death rate was p r o b a b i l i s t i c . I t was written i n IBM 360 assembler so that each colony was handled as a separate sub-model. The r e s u l t s from that model are described i n a l a t t e r chapter. 18 l i . . THE SET - THE TOOLS FOB THIS PARTICULAR PROBLEM Once an application i s decided on a dialogue should begin between the computer scientist and the biologist who has studied that particular system. There i s a logical set of questions that must be answered before the model can be constructed., In this chapter, these questions w i l l be taken in order and the ramifications of the answers pointed out in addition to explaining what the answers were i n this particular study* The processes of questioning and dialogue i s obviously an iterative one., Following an i n i t i a l set of discussions with W. Thompson, during which we touched on a l l these issues, I proceeded to read publications on the Western Tent Caterpillars which brought up further points needing c l a r i f i c a t i o n . These were discussed with W. Wellington. Finally, I had a long afternoon with the chief research person, P. Cameron.„• Although there were several iterations, the following questions remain the basis for discussion. 19 3.1 SPATIAL BOUNDARIES The f i r s t question to be answered i s : what are the s p a t i a l and time boundaries of the system? In some ways t h i s begs the l a t e r question of what are the independent variables, but as a rule i t turns out to be either s p a t i a l coordinates or time. Of course the system may be a subsystem of some larger system - say the universe - but reasonable boundaries can usually be drawn by looking at the environment of the e n t i t y or process under study.. For the c a t e r p i l l a r study, the s p a t i a l range was the Saanich peninsula and so s p a t i a l l y , the choice was rather easy.. Since the model .is interested i n o v e r a l l trends, the time i n t e r v a l has to be ,rather large, a factor several times larger than the basic generation time of one year. 3.2 COMPONENTS OF THE SYSTEM The second question i s answered by looking at the flow of energy or mass within the system and asking f o r the places where that energy or mass seems to pool* In the c a t e r p i l l a r model, we see the c a t e r p i l l a r s themselves, t h e i r food source, t h e i r predators, the weather, t h e i r tents, and the c a t e r p i l l a r colony locations as the major components. Another way to view components i s by the "process" which they are a part of. Bhat do c a t e r p i l l a r s do? They are l a i d as eggs, hatch, t r a v e l to and from egg mass to food source, eat, b u i l d a tent, r e s t , are eaten, expose themselves to weather, pupate, disperse and l a y eggs* That, of course, i s a very s i m p l i f i e d statement of t h e i r a c t i v i t i e s . In actual f a c t , even t h e i r resting . i s a 20 complicated process, where the temperature, food intake, individual caterpillar characteristics, and position relative to each other, determine the length of the rest and consequently when the next feeding begins. , In Chapter 4, I discuss how this l i s t of components was pared down to a workable size of major elements in the system.,,. 3.3 STATE AND DRIVING VARIABLES From the l i s t of system components, the third question regarding state variables arises.. A state variable i s one that uniquely determines one portion of the system*s status for any given value of the independent variable*. Hence a l l of the state variables together uniquely determine the particular state of the system as a whole at that particular instant., These state variables (or their derivatives) must be a function of some other variables in the system. In that sense, they . are dependent variables. The only variable mentioned above for the caterpillar model that was not a state variable was weather. . The next question to be considered i s that of driving variables. A driving variable i s one which i s not a function of any other system variable except the independent variable. In this model, weather i s the only driving variable. As I discuss in Chapter 4, weather i s thought to be the strongest single influence on caterpillar activity and population size., Some systems have no driving variables - that i s , each variable i s a function of at least one variable besides time., 21 3.4 THE INDEPENDENT VARIABLE In determining the independent variable, i t i s necessary to determine over what range the model's behavior i s to be studied. , As mentioned e a r l i e r , t h i s i s usually some period of time, but i t can take other forms, such as a measurement of distance. , 3.5 CONTINUOUS OR DISCRETE Probably the guestion with the most ramifications i s the one of continuity of the system being modeled. , A continuous system i s one i n which the change i n the system i s on going at every instant, where i n a discrete system, the change happens instantaneously at p a r t i c u l a r points i n time.,• I t has been said that the continuous i s but the l i m i t of the discrete and that the discrete i s but the time s l i c e picture of the continuous. That i s , a system may involve the appearance and disappearance of i n d i v i d u a l s in the system, but as the number of i n d i v i d u a l s increases, the e f f e c t of one i n d i v i d u a l a r r i v i n g or leaving the system approximates more and more a smooth continuous curve. , Inversely, i f one takes a smooth curve describing the population of a system over time, taking a sample of the system every change i n time T where T becomes large, becomes a bar graph of i n d i v i d u a l s entering and leaving the system. , A l l that i s to say that the choice of continuous against discrete modeling depends more upon the magnitudes of the numbers involved, rather than the types of e n t i t i e s being modeled. , A rather controversial h e u r i s t i c can be employed to make 22 t h i s decision of continuity. With a few exceptions, applications i n the b i o l o g i c a l and physical sciences would use continuous system . simulations while those i n commerce, transportation, and s o c i a l sciences would probably - use discrete system simulators. Additional i n s i g h t i n t h i s matter can be gathered by comparing the advantages of continuous and discrete simulations f o r the p a r t i c u l a r application i n question. Discrete systems are e s p e c i a l l y h e l p f u l i n developing s t a t i s t i c a l information about the system: average delay, maximum and minimum delay, average process time etc., Discrete systems also deal with the well defined, complex, interactions between i n d i v i d u a l s i n a system where "events", " a c t i v i t i e s " , or "processes" take place. The important f e e l i n g one should have about "events", " a c t i v i t i e s " , and "processes" i s that they should be a well defined function between in d i v i d u a l s i n the system. In such a case, discrete simulation i s the best approach.., Continuous systems, on the other hand, deal with flow rates between segments of the system, and are not so concerned with s t a t i s t i c s of the system, as they are with o v e r a l l numbers i n a given component. Continuous systems can deal with a large number of i n d i v i d u a l s , because the population i s treated as a single number - the t o t a l population. Discrete systems can not deal with large populations, because they must handle each i n d i v i d u a l separately and for a large number, t h i s reguires too much time., For the c a t e r p i l l a r model, the problem involves a large 23 number of individuals (approx., 1,000,000) and the system responses that are of interest are overall numbers of individuals. J These two conditions caused me to choose a continuous model. As i t turned out, one of the major features of the real world system was a periodic smooth curve that described the population over a long time span, so one naturally began to think in terms of continuous simulators. The primary result of deciding on the continuity of the system, i s the range - of languages available for the simulation., If a discrete simulation was chosen, one could either use a general algorithmic language li k e ALGOL or SIMULA or a specialized discrete simulation language like GPSS. For the continuous system, the specialized language would be something like CSHP or BEDSOCS. , Although these languages are more powerful than their predecessors, they are not available everywhere and are not understood by the majority of the s c i e n t i f i c community* As I wanted to be able to run this program in any computing environment and allow others to use and modify i t , I decided to use FORTRAN, which i s far and above the most common sc i e n t i f i c programming language, to simulate the continuous caterpillar system* , 24 3.6 STOCHASTIC ELEMENTS The f i n a l question relates to the stochastic nature of the system or the element of chance involved. For example, the rate of a r r i v a l of t r avelers at a border crossing i s not constant, but random. The rate i s therefore stochastic i n nature and has a measure of uncertainty i n i t . This guestion i s s i m i l a r to the question of continuity i n that the decisions are a r b i t r a r y rather than cl e a r cut. Certainly i n the r e a l world, almost everything has an element of chance i n i t . Which deer a mountain l i o n may k i l l i s , a f t e r the set of possible victims i s established, a matter of which one i t "comes upon" when i t i s hungry* The saving grace i s that when large numbers are involved, one can sometimes ignore the stochastic elements*, I f most deer behave somewhat s i m i l a r l y , i t doesn't matter which one of a large number of deer the mountain l i o n eats. The r e s u l t i s the same - the deer population decreases by one. Hence one could say one of two things: the p r o b a b i l i t y of one deer dying on a given day i s . 0 1 ; the rate of mortality i s approximately . 0 1 deer/day. The deterministic statement i s usually f a s t e r i n simulation, and so i f the difference between the two i s small (as i t was f o r the c a t e r p i l l a r s ) , the deterministic model may be used. . The f i r s t (and. more detailed) model of tent c a t e r p i l l a r s [ 7 ] by W. Thompson et a l had many stochastic elements. One a d d i t i o n a l clue i n deciding on the stochastic c h a r a c t e r i s t i c s i s to see i f the difference introduced by the chance element i s magnified by some gross change i n the state 25 at that point. Consider for a moment t h i s problem when the females disperse at the end of the year. A l l year long the c a t e r p i l l a r s are r e l a t i v e l y stationary and reducing i n numbers from several hundred per colony to a fev per colony. Those c a t e r p i l l a r s pupate and f l y i n some ar b i t r a r y (?) d i r e c t i o n , some a r b i t r a r y (?) distance, and deposit t h e i r eggs - a l l 200 of them. , Hence 200 c a t e r p i l l a r s may be placed i n any gr i d sguare. That large number of c a t e r p i l l a r s "magnifies 1' the e f f e c t of chance i n the migration of the moths to the point that I could no longer use a purely deterministic model. The migration d i r e c t i o n and distance had to be stochastic. A f i n a l personal bias - not e n t i r e l y unfounded - i s added here., As I have noted above, the r e a l world i s almost always discrete and stochastic., Consider the atom and random electrons which are the building blocks of matter to convince yourself of that f a c t . Nevertheless, i t appears to me, a reasonable h e u r i s t i c i n modeling i s that to s t r i v e f o r s i m p l i c i t y , continuity, and deterministic modeling. A optimistic engineer has been known to say "Nature i s , on the whole, deterministic, continuous, and twice d i f ferentiable.'? 26 Us. THE SCRIPT - IN RETROSPECT I emphasize the t i t l e of this chapter. As i s pointed out by Wellington 1s paper on approaches in population dynamics [12] , the work published on a given subject area i s almost always polished to the point that one reading i t might assume that things naturally f e l l into place a l l the way along during the research and that the course of action was clear from the beginning of the investigation. That of course i s never the case. If necessity i s the mother of invention, then surely serendipity i s the father. When I f i r s t decided upon caterpillars as my . topic, I was going to study a detailed mortality model of the system, but as I pursued the various factors involved with mortality (such as random attack, starvation, non-random . attack, desiccation and v i r a l infection), I decided that the area that needed studying was a more general mortality model which had only a few terms dependent on the major state variables. So the work that resulted i s seen here i n retrospect not as a preset goal. In any application of modeling I would see a dynamic direction of effort that would change as new factors and understanding come to bear. One would begin to assemble the answers to the questions in the previous chapter, and as the model takes shape, one would see questions that require further study. Perhaps the whole model would be scrapped in favor of studying an interesting subsystem. 2 7 4.1 DETERMINATION OP THE MAJOR STATE VARIABLES This chapter's text naturally f a l l s into three t o p i c s : i ) the determination of rate equations f o r the major state variables, i i ) the detailed description of the c a t e r p i l l a r model, and the determination of the parameters i n those rate equations. . Let me review i n d e t a i l the sequence of events that lead me to choose d i f f e r e n t i a l equations to model the relationships i n the tent c a t e r p i l l a r community. There i s some data that could be used to construct some complex and well defined r e l a t i o n s h i p s between i n d i v i d u a l c a t e r p i l l a r s that e f f e c t t h e i r a c t i v i t y , eating habits and loss to predators. , As an example of t h i s , there was a study done on four colonies of c a t e r p i l l a r s . , Two of the colonies had a high percentage of active i n d i v i d u a l s and the other two had a l o w percentage of active i n d i v i d u a l s . , One of each type was placed i n a harsh weather s i t u a t i o n and one of each type was placed i n a mild weather s i t u a t i o n . Detailed records were kept on the types of m o r t a l i t i e s ; predation by spiders, desiccation, etc.. From that kind of data one might be lead to construct a very detailed mortality model, and since the data on these relationships i s f a i r l y r i c h , one might even pursue a model for a i n d i v i d u a l colony. . This could even be a discrete model, i f only a single colony was studied* Two f a c t o r s kept me from t h i s approach.. Primarily, I was interested i n the c a t e r p i l l a r population as a whole (which as pointed out e a r l i e r , was f a r too large f o r a discrete simulation). That was more of a personal preference, but when coupled with the fact that 28 weather was f e l t to be the overriding factor i n c a t e r p i l l a r population dynamics, the scales were tipped i n the d i r e c t i o n a macro, continuous model. In my discussions with Thompson, Cameron and Wellington, weather was alluded to again and again as the major c o n t r o l l i n g element. I t i s important to see that the b i o l o g i s t was the one to make t h i s judgement. The computer s c i e n t i s t i s required merely to confirm the supporting arguments fo r t h i s judgement. In the case of c a t e r p i l l a r s , the supporting evidence was more along the negative l i n e s as to which things did not e f f e c t the population. In most b i o l o g i c a l . systems, food supplies and predator populations regulate the prey population, but i n t h i s system i t was not the case. , Even i n the most populated years, the food supply was almost unlimited. In a very few cases did a c a t e r p i l l a r colony t o t a l l y d e f o l i a t e i t s host t r e e . There was some starvation, of course, but that appeared to be from the i n a b i l i t y to forage from food. There were also a few c a t e r p i l l a r predators.. Wasps, spiders, ants and several in s e c t parasites a l l prey on c a t e r p i l l a r s , but usually they did not have a s i g n i f i c a n t e f f e c t . , Besides that, the predator and parasite populations were nearly independent of the prey populations because they had alternate hosts i n years of low prey populations. The general f e e l i n g , then was that weather played the biggest role i n c o n t r o l l i n g c a t e r p i l l a r populations., Weather i s a d i f f i c u l t thing to analyze quantitatively., 29 Most of i t s e f f e c t s are i n d i r e c t . Bad weather may keep c a t e r p i l l a r s i n t h e i r tents away from food and thereby k i l l them by starvation. Good weather may allow the c a t e r p i l l a r s to b u ild new tents, leaving the weaker, virus infected members of the colony behind i n the previous tents and i s o l a t i n g what could have been an epidemic. Again, the e f f e c t i s i n d i r e c t . Because these e f f e c t s are i n d i r e c t , one can only model them i n general, and not by a d i r e c t , measurable r e l a t i o n s h i p . This i s where d i f f e r e n t i a l eguations become quite helpful* Suppose a simple rel a t i o n s h i p between weather and the mortality could be expressed i n d i f f e r e n t i a l eguations that would f i t the data for actual c a t e r p i l l a r systems. Even though the d i f f e r e n t i a l equations would be a r t i f i c i a l , i t could express these i n d i r e c t relationships and would apply . the e f f e c t uniformly i n the system over time. , As noted e a r l i e r , the o v e r a l l system responses seem to be continuous. , So X had answered two basic questions: the type of equations and the dri v i n q variable. Next, I looked back over the early work of Wellington., Surely . his biggest contribution to biology was that i n d i v i d u a l s behavior, or rather types of i n d i v i d u a l s and t h e i r behavior, a f f e c t the o v e r a l l c h a r a c t e r i s t i c s of the ecosystem. The active i n d i v i d u a l s i n the c a t e r p i l l a r colony lead the re s t of the colony to new foraging areas* .. They are the f i r s t to terminate the r e s t periods of the colony. They also influence the number and shape of the tents constructed by the colony. Their dispersal distances as moths are much longer. , The sluggish members of the colony retard i t s a c t i v i t y , but they 30 do reinforce i t s tent construction., What i s the simplest statement that could be used to model these differences? One p o s s i b i l i t y would be to divide the c a t e r p i l l a r s into groups by type., In the studies on the Saanich Peninsula, t h i s i s exactly what was done. , The four a c t i v i t y groups decided upon were r e a l l y . very a r t i f i c i a l . In r e a l i t y , a spectrum of a c t i v i t y existed, ranging from those that were not even able to break out of t h e i r eggs, to those that were very independent and constantly active. S t i l l , an important point can be made here. When faced with a population that i s guite d i v e r s i f i e d i n i t s behavior, one must t r y and divide i t into prototypes and then model these groups of in d i v i d u a l s as separate populations, . so that the differences i n in d i v i d u a l s can be incorporated., This model required four rate equations, one for each type of i n d i v i d u a l * The four populations corresponding to these four types became the state variables of the system*. Once the differences were established, the int e r a c t i o n between them was chosen to be a constant times the percentage the number of most active i n d i v i d u a l s were of the whole* , The one remaining factor to be incorporated i n the state equations was the independent variable, time. , I t turns out that the mortality rate i s not constant. C a t e r p i l l a r s shed t h e i r skin, and t h e i r age (instar) i s measured by the number of times that has occurred. The mortality rate i s quite high in early i n s t a r s but as the c a t e r p i l l a r grows, i t becomes l e s s susceptible to many predators and diseases, such as ants which 31 can only attack young c a t e r p i l l a r s . In the l a s t i n s t a r , the c a t e r p i l l a r becomes vulnerable to f l y parasitism. , So the natural thing i s to add constant m u l t i p l i e r s and time to the state eguations. , I t was then possible to write the simplest d i f f e r e n t i a l equations that express the four state variables i n terms of themselves, weather ( I ) , and time (T). P(i) i s the population of the i t h type of c a t e r p i l l a r . The vectors A,B,C, and D are the parameter values to be determined l a t e r . This gives the following set of equations. dP(i)/dT = (A(i) • B(i)*T + C(i)*W + D (i) *PCT) *P (i) where i goes from 1 to 1 PCT i s p(1)/S and S. i s the sum of the P(i) These four equations became the hub of my model. Routine INTGBT takes the current values of the state variables and time and c a l c u l a t e s the derivatives d i r e c t l y from these equations. That routine could be c a l l e d by any integration scheme to c a l c u l a t e the state of the system at a new time T» = T. + DT. , In the f i r s t stages of development, I used Eulers integration method to numerically integrate a day at a time. See Appendix I for a description of the Euler algorithm., v 32 4.2 DESCRIPTION OF THE MODEL He now return to the spatial considerations discussed i n chapter III., I wanted to study the spatial distributions of caterpillars on the Saanich Peninsula. The approach used was very similar to the one used to break . down the caterpillar population into prototypes. Suppose the peninsula was to be broken down into grid squares, and each grid square was assumed to be homogeneous in micro climate, tree height, and general caterpillar occupancy., Each grid square would then have a different set of four populations of caterpillars, but the same four differential equations, except that the weather in each grid square would be altered by a function of the micro climate before using i t in the state differential equations. , Routine IHTGRT, then, .takes the current values of the four caterpillar populations and the local weather for one grid square and returns the derivatives which are integrated ahead one day at a time for a l l the grid squares in routine STEP. , The grid squares with nonzero populations are kept in a double linked l i s t so that each square need not be searched each day to see i f i t must be integrated. The main program i s quite straight forward. - It call s READ to input a l l the i n i t i a l population data.. It call s INITAL to i n i t i a l i z e a l l the variables, and then calls STEP for as many days as i s required to integrate the system to the end of the year. DSPRS i s called to change the remaining caterpillars to moths, disperse them randomly, and lay new egg masses. ,• A new population distribution i s then printed out by OUTPUT and the 33 process begins again u n t i l the number of years of simulation required are produced., DSPSS and OUTPUT actually require * further explanation., The problem i n DSPHS i s that by t h i s time the c a t e r p i l l a r s have pupated and turned into moths, so they need to f l y i n some random d i r e c t i o n and distance and then deposit a quantity of eggs. , One soluti o n i s to pick a random number between 0.0 and 360.0, ca l c u l a t e the s i n and cos of that number, and go i n that d i r e c t i o n a random distance whose range of values i s determined by the type of i n d i v i d u a l (note the in d i v i d u a l s e f f e c t again). Even the number of eggs l a i d and the types of ind i v i d u a l s r e s u l t i n g from them are dependent upon the prototype of the mother.. That solution i s s u f f i c i e n t , but i t can be speeded up by noting that the maximum distance covered by a moth under any circumstances i s twenty grid sguares. Hence, with a f i n i t e number of equal arched rays, each gri d sguare could be v i s i t e d . Picture a bic y c l e wheel on a checker board* , I f enough spokes are added, each square would be touched by. at least one spoke* . It turns out, that the magic number of "spokes" i t takes to cover the twenty sguares i s twenty eight. Thus, only that many sines and cosines allow the moths to v i s i t any square within i t s range of f l i g h t . This f a c t saves the ca l c u l a t i o n of a s i n and cos f o r each moth (remember, there are sometimes several hundred thousand of them!) each year! I t was decided that i f a moth lands i n the ocean or on an uninhabited g r i d square, i t i s assumed to die. OUTPUT allows the user to output the yearly c a t e r p i l l a r 34 d i s t r i b u t i o n , either as a printer two dimensional plot with each g r i d position given a character from 1 to Z depending upon i t s t o t a l population, or as a three dimensional graphic r e l i e f map with the " a l t i t u d e " of the r e l i e f corresponding to the c a t e r p i l l a r population. I f the r e l i e f map i s displayed continuously as the population changes from year to year, one begins to get a vi s u a l f e e l i n g f o r the behavior of the system over a long period of time. 4.3 DETERMINATION OF THE PARAMETER VALUES The next step of model development was the most d i f f i c u l t one: the determination of reasonable values f o r the parameters A, B, C, D, and E i n the d i f f e r e n t i a l equations describing the population rate of change., As pointed out e a r l i e r , these parameters could not be determined from d i r e c t r e l a t i o n s h i p s in the f i e l d data. , Heather and the e f f e c t s each type of c a t e r p i l l a r has on the other types are i n d i r e c t e f f e c t s . I decided to take a couple of average colonies i n average weather and harsh weather and determine the parameters that would r e s u l t i n the f i n a l colony populations being as close to the observed values as possible., That meant that I had a non-l i n e a r system of d i f f e r e n t i a l equations to be solved as a two point boundary value problem. Given an i n i t i a l set of conditions, and the d i f f e r e n t i a l equations as a function of the parameters A,B,C, and D, the parameters are varied so that the r e s u l t i n g f i n a l population meets a given terminal condition. The number of degrees of freedom was equal to the number of parameters (twelve)., The number of i n i t i a l and 3 5 terminal conditions was the number of state variables, times the number of test colonies (eight). Actually, I even experimented with three test colonies, but that meant the solution (if i t existed) was quite d i f f i c u l t to find., After some effort, I reduced the number of test colonies to two. At f i r s t , I took the simple minded approach of ••twiddling1' each of the parameter values until the c r i t e r i a were nearly satisfied. , It became like stuffing a f u l l inner tube into a bicycle t i r e . As soon as I poked i t in one place, i t popped out in another., I needed something more powerful* The second try was with a Newton-Raphson iteration technique. Given a error vector (E) which i s how far away each terminal condition was from a given i n i t i a l guess of the parameters, and a sensitivity matrix (S) which i s the partial derivatives of the errors with respect to the parameters, one can compute a correction vector 0 = -S 1(SS 1)~*E. This can be used to compute a new set of parameters, a new error vector, and a new sensitivity matrix., The process is. used iteratively until the errors are within some tolerance of zero* To calculate E for any given values of the parameters, one needs numerically integrate the state differential equations for one year. , To calculate S for given values of the parameters (P)i one takes the nominal values P(0) and perturbs them one at a time and calculates the corresponding E vector. / S i s then the values (E(0) - E(j))/(P(0) - P(i)) where j varies over the number of errors and i varies over the number of parameters. The procedure begins by calculating E(0) from the i n i t i a l 3 6 guess on the parameters P (0). That requires integrating the d i f f e r e n t i a l . equations forward to the terminal time and subtracting the terminal population values from the desired ones., Next, one at a time, the parameters are increased by a small constant value and the d i f f e r e n t i a l equations are integrated to the terminal time where the new error values E(j) are calculated. As each parameter i s perturbed and the corresponding E (j). i s calculated, a new row i s added to the se n s i t i v i t y , matrix S. Ihen the process i s f i n i s h e d , the vector U i s calculated and new parameters are calculated from P,=P+0., A new error vector E(0) i s then calculated by integrating the d i f f e r e n t i a l equations again. I f i t i s within some tolerance of zero, the process i s completed and P* i s printed out., Otherwise, P* becomes P and the process begins again. This Newton-Raphson approach appeared to work quite well, u n t i l I examined the intermediate values of the derivatives of the state variables. To my horror, they were posit i v e at some points i n . time. „ Mathematically, that was quite reasonable, but b i o l o g i c a l l y , c a t e r p i l l a r s do not reproduce or spontaneously . appear during the year. The population must always decrease from the time the eggs are l a i d u n t i l the time when they pupate and disperse at the end of the year. Hy f i r s t modification was to force the derivatives to be negative by including the negative of t h e i r absolute value. That turned out to be too a r t i f i c i a l , and the Newton-Raphson scheme began to diverge. I r e a l l y needed some way of 37 constraining the d i r e c t i o n of the correction vector and some intermediate values of the derivatives. I f I had wanted equality constraints, I could have just added them to the l i s t of eight I was already targeting f o r , but I needed inequality constraints., I didn't want to constrain the derivatives unless they started to go pos i t i v e ^ and I didn't want to constrain the correction vector unless the correction d i r e c t i o n would v i o l a t e known b i o l o g i c a l f a c t s . These inequality constraints forced me to use a program developed e a r l i e r with some engineers i n trajectory shaping£2]. Now I had the t o o l s necessary to attack the problem. after modifying t h i s program extensively and playing with the constraints so that the derivatives behaved properly, I came up with the following values f o r the parameters. A(i) . = .009960 .120000 .30700 .009350 B(i) = -.000685 -.000126 -.00055 -.000542 C(i) = -4000151 -.003850 -.01000 -.000289 D(i) = -*001000 -.001000 -.00100 -.001000 The reader i s referred to Appendix II for the source code f o r the Newton-Raphson routines and Appendix III f o r the inequality constraint programs. , This Figure I and Figure I I show . the time history of both the state variables and t h e i r derivatives. , FIGURE I - STATE VARIABLES IS AVERAGE HEATHER FIGURE I I - DERIVATIVES OF STATE VARIABLES IN AVERAGE HEATHER no 5 . VALIDATION One of the hardest things to accept as a computer s c i e n t i s t was the statement that when a model produces r e s u l t s that " f i t " known observed f a c t s , that i t hasn't proved the correctness of the model. , In f a c t , a model can never be shown to be an exact mapping from the o r i g i n a l system, no matter how s i m i l a r i t may behave. No matter how consistently a model seems to duplicate known f a c t s , i t has not proved i t s correctness. The best one can do i s make a good case f o r the l i k e l i h o o d of the model being a reasonable one. To an absolute and objective computer s c i e n t i s t , that i s humbling r e a l i z a t i o n ! ! . One can, however, show a model to be incomplete or erroneous by a counter example where the r e a l system behaves d i f f e r e n t l y than the model., Even though i t i s not possible to prove the correctness of a model, there are some ways to validate i t or show the l i k e l i h o o d that i t w i l l behave i n most cases as the r e a l system would behave. These methods of vali d a t i o n must a l l be non-circular; that i s , the data that they use f o r comparisons must not have been used i n constructing the model.. I f one uses the temperature history of the Fraser r i v e r from 1950 to 1975 i n constructing a salmon modeling program, he would not be able to validate his model by using the comparison of how the model . behaved to how the r e a l system behaved i n 1960, since that was part of,the data he used to bu i l d or modify h i s model i n the f i r s t place. Ask someone a question. Now give them the answer., Now ask them the question again, and i t i s 41 surprising how many times that person w i l l have the answer! The data used i n val i d a t i n g a model must be as unrelated as possible to the data used to construct the model i n the f i r s t place., Validation i s usually done by comparing the state of the model a f t e r some long change i n the independent variable with the r e a l system. ,. As noted before, t h i s i n t e r v a l should be di f f e r e n t than any i n t e r v a l used i n " f i n e tuning" the model. Additional v a l i d a t i o n i s possible by comparing the model's reaction to changes i n the driv i n g variables with the reaction of the r e a l system. , This type of va l i d a t i o n i s sometimes c a l l e d s e n s i t i v i t y analysis. It centers upon the o v e r a l l trends i n and the r e s i l i e n c e of the system and i t s model. Resilience i s the a b i l i t y of the system ;to return to a steady state a f t e r i t i s perturbed by some external force. In the c a t e r p i l l a r model, both types of vali d a t i o n were used. Figure III shows the model's c a t e r p i l l a r population, the Thompson et a l model's r e s u l t , and that of the r e a l population from 1959 to 1963. , 42 FIGUBE I I I - COMPABISON OF HEAL AND SIMULATED POPULATIONS These years of data were not used i n building the model. The two were hard to compare, because the f i e l d study count was i n number of colonies and the model's output were i n actual c a t e r p i l l a r s . The number of c a t e r p i l l a r s per colony varies from year to year. Generally the number of c a t e r p i l l a r s per colony increases as the weather improves from 40 to 90. The number of colonies might be wrong by a factor of two. Considering the s i m p l i f i e d nature of the state equations, these were remarkably close to the r e a l values. 43 The second area of validation i n the c a t e r p i l l a r model involved the model's- reaction to weather. During the c a l c u l a t i o n of the system parameters discussed i n chapter four, only two constant weather values had been used. Figure IV shows the r e a l weather and the response of the model to i t from 1959 to 1964. 44 C 3 . IT. r-' o in" in ra 1959.0 1960.3 1961.5 1962.8 1964.0 FIGURE IV - WEATHER AMD THE SIMULATED CATERPILLAR POPULATION Notice that when the weather improved from 1959 to 1963, the population increased, and i n 1964 when the weather was harsh, the population decreased. A cl e a r e r picture of what happens s p a t i a l l y can be found by p l o t t i n g the c a t e r p i l l a r population on a two dimensional g r i d . . Figure V i s a plot of the population on the Saanich Peninsula i n 1960 when the weather was quite mild. FIGURE V - CATERPILLAR DISTRIBUTION IN HILD WEATHER Figure VI i s a plot of the population i n 1964 when the weather was quite harsh* FIGURE VI - CATERPILLAR DISTRIBUTION IN HARSH HEATHER Note that the milder weather allowed the c a t e r p i l l a r s to i n f e c t areas with r e l a t i v e l y h o s t i l e l o c a l climate, but i n harsh weather, only the areas with the best l o c a l climate were occupied. The refuge areas into which the c a t e r p i l l a r s receded i n years of harsh weather agrees with f i e l d observations. Figures VII and VIII demonstrate another method of p l o t t i n g the mild and harsh weather d i s t r i b u t i o n s . 47 FIGURE ¥11 - CATERPILLAR DISTRIBUTION IH MILD WEATHER 4 8 FIGURE YIII - CATERPILLAR DISTRIBUTION IN HARSH WEATHER 49 6 . CONCLUSIONS One of the most encouraging things that happened, was that with relatively simple state equations expressed as ordinary differential equations, there was reasonable agreement with the real ecosystem. , The simplified nature of the model i s quite appealing in that i t makes i t s functioning quite straight forward and allows for the modification of the state differential eguations. , Furthermore, the graphics display of the caterpillar population on a two dimensional grid (which produced figures V and VI ) gave a visual understanding of what i s taking place. „ Currently, work i s being done on taking a movie of such displays over a long period of time, and should give added appreciation for the functioning of refuge areas and the response of the model to varied weather conditions. Areas of future study would undoubtably include extending the simple differ e n t i a l equations to include additional terms, such as cross product terms between weather, time, and percent of type I caterpillars. Perhaps an extension such as E(i)*W*T + F(i)*H*PCT would be a good place to start.. One would tend to expect better agreement with the real ecosystem as more detail i s added to the differential equations in the model., Of course new terms in the differential equations would necessitate a new set of parameter values. Experience has shown that the determination of those values should not be limited to examining a single colony in two weather conditions, but should be-done using a whole grid square of 50 colonies at one time. , Another area that would be i n t e r e s t i n g to address, involves long periods of hypothetical weather. Using the s i x years of known weather to create r e a l i s t i c long term weather, i t would be productive to study extended periods of harsh and mild weather. The r e s u l t s would show how responsive the model i s o v e r a l l to weather and might be used to predict the conditions under which extinction or saturation might occur.. The v i s u a l plots produce by the model hold the greatest prospect for future use. The movie included with t h i s t h e s i s demonstrates very c l e a r l y the important part v i s u a l images can play i n understanding the functioning of a simulated system. That discovery i s not only relevant to bio-simulation, but i s a powerful technique to be applied to modeling i n general. 51 7. , B I B L I O G R A P H Y [I] Bucher, G.E., "Winter Rearing of Tent Caterpillars, Halacosoma ", Canadian Entomologist, Vol. 91, 1959, pp. 411-416., [2] Cornick, D.E. and R.T. Steinhoff, "Uphill, a Constrained Optimization Program", (unpublished), 1973. £3] Hay, R.H., Stability and complexity i t Model Ecosystems, Princeton University Press, Princeton,1973. [4] Sammets, J-E«# Programming Languages; . History and Fundamentalsff Prentice-Hall, Englewood C l i f f s , 1969, pp. ,627-631, 650-660. [5] "The SCi Continuous System Simulation Language (CSSL)", SimulationVol. .9. Ho..6, 1967, pp.,281-303. , £6] Sullivan, CR.and H.G* Wellington, "The Light Reactions of Tent Caterpillars Malacosoma Pluyiale (Dyar)", Canadian Entomologist. Vol. 85, 1953, pp..297-310. , [7] Wellington, H.G., P..Cameron, W.A. Thompson, I.,Vertinsky, and A. Landsberg, "A Stochastic Model for Assessing the Effects of External and Internal Heterogeneity on an Insect Population,, Resource Population Ecology. Vol. ,17, 1975, pp. 58-85. [8] Wellington, W.G., "Individual Differences as a Factor in Population Dynamics: the Development of a problem", Canadian Journal of Zoology a Vol. 35, 1957, pp. 293-323. [9] Wellington, W.G., "Individual Differences in Larvae and Egg Masses of the Western Tent Caterpillar'?, Canadian Department of Agriculture Forest Biology Division Bi-monthly Progress Report^ Vol. ,-15, No. 6, 1959. [10] Wellington, W.G., "Qualitative Changes in Natural Populations During Changes in Abundance", Canadian Journal of Zoology, Vol. ,38, 1960, pp. 289-314. [II] Wellington, W.G., "Qualitative Changes in Populations During in Unstable Environments", Canadian Entomologist. Vol. ,96, 1964, pp..436-451. [12] Wellington, W.G., "An Approach to a Problem in 52 Population Dynamics", Forest Entomology and Pathology i Vol. 1, 1965, pp* , 175-185. APPPENDIX I - NUMERICAL METHODS EULERS METHOD 53 The general eguations 1) D(i) = dP(i)/dt = F(i) where F i i s a function 2) P« (i) = P(i) • D(i)*(T-T0) of P and T The general process i) calculate vector D from eguation 1 i i ) calculate vector P* from eguation 2 i i i ) P* now becomes P and TO becomes T APPPENDIX II - NEWTON-RAPHSON ROUTINES MAIN PROGRAM 54 1 IMPLICIT REAL*8 (A-H,0-Z) 2 DIMENSION X(12), G(10), GX(10,12), THIX(3), DX 1 (3) , DX2 (12) 3 DATA THI /O.ODOO/ 4 N = 0 5 NF = 0 6 NC = 10 7 FACT2 = .5 8 READ (1,5) X 9 WRITE (6,3) X 10 CALL STATE (10,12,NC,X,G,GX,GNEW) 11 C WRITE (6,3) FACT2 12 WRITE (6,2) GNEW 13 WRITE (6,3) (G (I) , 1 = 1,NC) 14 DO 9 J = 1, 12 15 9 WRITE (6,3) (GX (I, J) ,1=1,NC) 16 10 CONTINUE 17 CALL DX (NC,12,10,G,GX,THIX,DX1,DX2,2) 18 WRITE (6,3) DX2 19 READ (5,5) FACT2 20 DO 101 I = 1, 12 21 X(I) = X(I) - DX2(I)*FACT2 22 101 CONTINUE 23 C WRITE (6,3) FACT2 24 WRITE (6,3) X 25 GOLD - GNEW 26 CALL STATE (10,12,NC,X,G,GX,GNEW) 27 N = N + 1 28 WRITE (6,2) GNEW 29 WRITE (6,3) (G (I) ,1=1,HC) 30 DO 102 J = 1, 12 31 102 WRITE (6,3) (GX (I,J) ,I=1,NC) 32 IF (GNEW .LT. GOLD) GO TO 11 . 33 FACT2 = FACT2/2 34 NF= 0 35 GO TO 20 36 11 CONTINUE 37 IF ( (DLOG10(GOLD)-DLOG10 (GNEW) ) .GT. 1.0) GO TO 12 38 NF = NF + 1 39 IF (NF .LT. 3) GO TO 20 40 FACT2 = FACT2*2 41 IF (FACT2 .GT.,1) FACT2 = 1 42 NF = 0 43 GO TO 20 44 12 CONTINUE 45 FACT2 = 1.0 46 NF = 0 47 20 CONTINUE 48 IF (GNEW .GT. 1.0D-10) GO TO 10 49 WRITE (6,1) (X(I) ,1=1,12) 50 STOP 51 1 FORMAT (1X, 6(D17.10,3X)) APPPENDIX II - NEHTON-BAPHSON ROOTIHES MAIN PROGRAM 55 52 2 FORMAT (10X,10H EBBOBMAG ,D14.7) 53 3 FOBMAT (1X,10(2X,D10.3)) 54 4 FOBMAT (12,12) 55 5 FOBMAT (D8*2) 56 END APPPENDIX II - NEWTON-RAPHSON ROUTINES SOBfiOOTINE STATE 56 1 SUBROUTINE STATE (L,H,N,XO,G0#GX,GMAG) 2 C 3 C CALCULATES THE STATE OF THE SYSTEM 4 C ' 5 C GO IS THE CURRENT ERROR 6 C GX IS THE CURRENT SENSITIVITY MATRIX 7 C 8 IMPLICIT REAL*8 (A-H,0-Z) 9 DIMENSION X(20), XO (H) , G(20), GO (N) , GX.(L, M) , GSUM(20) 10 DATA NTIMES /-1/ 11 DO 10 I = 1, II 12 GSDM (I) = 0 13 10 CONTINUE 14 DO 20 I = 1, M 15 X(I) = X0(I) 16 20 CONTINUE 17 CALL ERROR(M,N,X0,G0) 18 GMAG = 0 19 DO 30 J = 1, N 20 GMAG = GMAG + DABS (GO (J) ) 21 30 CONTINUE 22 WRITE (6,1) (GO (I) ,I=1,N) 23 NTIMES = NTIMES + 1 24 C IF (M0D(NTIMES,4) .NE. ,0) RETURN 25 DO 100 I = 1, B 26 X(I) = X0(I) • 1.QD-2*X0(I) 27 CALL ERROR (M,N,X,G) 28 DO 40 J = 1, N 29 GX(J,I) = (G(J) - G0(J))/(X(I)-X0(I)) 30 GSDM(J) = GSUM(J) + DABS (GX (J, I) ) 31 40 CONTINUE 32 X(I) = XO(I) 33 100 CONTINUE 34 DO 120 J = 1 r N 35 GSUM(J) = GSUM(J)/N 36 IF (GSUM(J) .EQ. ,0) GSUH (J) = 1.0D0 37 GO(J) = G0(J)/GSUM(J) 38 120 CONTINUE 39 DO 160 J = 1, N 40 DO 140 I = 1, M 41 GX(J #I) = GX(J#I)/GSUM(J) 42 140 CONTINUE 43 160 CONTINUE 44 RETURN 45 1 FORMAT (1X,10(2X,D10.3) ) 46 END ( APPPENDIX II - NEWTON-RAPHSON ROUTINES SUBROUTINE ERROR 57 1 SUBROUTINE ERROR (H,N,X,G) 2 C 3 C CALCULATES THE ERROR BI INTEGRATING THE 4 C EQUATIONS TO THE FINAL TIME 5 C 6 IMPLICIT REAL*8 (A-H,0-Z) 7 COMMON /COM/ C (12) 8 DIMENSION X(M). G (N) , P(13), P0(13), PD(13), Q{13) 9 DATA PO /O.ODO, 10 1 5.0D1,7.1D1,1.12D2,3.0D0, 11 2 5.0D1,7.1D1,1.12D2,3.0D0, 12 2 2.6D1,4.5D1,9.9D1,1.5D1/ 13 DO 10 I = 1, 13 14 P(I) = P0(I) 15 10 CONTINUE 16 DO 20 I = 1, M 17 C(I) = X(I) 18 20 CONTINUE 19 DO 40 I = 1, 80 20 CALL DRK(P,PD,Q,1.0D0,13,1) 21 DO 30 J = 2, 13 22 C IF (P(J) .LT. .1) P(J) = 0.OD0 23 30 CONTINUE 24 C WRITE (6#123) P 25 123 FORMAT (1X# 13010.-3) 26 40 CONTINUE 27 G (1) = P (2) - 3.0D0 28 G(2) ,= P(3) - 3.6D1 29 G (3) = P(4) - 1.1D1 30 G(4) = P(5) •- 7.OD0 31 . G(5) = P(6) - 1.5D0 32 G(6) = P(7) - 0.0D0 33 G(7) = P(8) - O.ODO 34 G (8) = P (9) - O.ODO 35 G(9) = 6.0D0-P(10)-P(11)-P(12) 36 G(10) = P(13) - O.ODO 37 RETURN 38 END APPPENDIX I I - NEHTON-RAPHSON ROUTINES SUBBOUTINE AUXRK 58 1 SUBROUTINE AUXRK(X,F) 2 C 3 C CALCULATES THE DERRIVATIVES (F) OF THE 4 C STATE VARIABLES (X) 5 C 6 I H P L I C I T REAL*8 (A-H,0-Z) 7 COMMON /COM/ C (12) 8 DIMENSION X (1) ,F (1) ,H (3) i P C T (3) 9 DATA 8 /2.5D1 #7.5D1,5.46D2/ 10 DO 20 I = 1, 3 11 14 = 4 * ( I - 1 ) + 1 12 T = O.ODO 13 DO 10 J = 1, 4 14 K = 14 • J 15 T = T + DABS(X(K)) 16 10 CONTINUE 17 P C T ( I ) = DABS (X (14 + 1) )/T - 2.0D-1 18 P C T ( I ) •= PCT (I) * 1 . 0D-1 19 I F (PCT(I) .GT. 0) PCT ( I ) = O.ODO 20 20 CONTINUE 21 DO 40 I = 2, 13 22 12 = I + 2 23 J = 12/4 24 K = 3 * ( I 2 - 4 * J ) 25 F ( I ) = (C(K+1) • C(K+2)*X(1) • C (K+3) * M (J) • PCT(J) ) *DABS (X (I) ) 26 I F ( F ( I ) .GT. 0) F ( I ) = 0 27 40 CONTINUE 28 RETURN 29 END APPPENDIX I I - NEHTON-RAPHSON ROUTINES SUBROUTINE DX 59 1 SUBROUTINE DX (NC,*NV,MC,G,GX,THIX,DX1,DX2, IS) 2 C 3 C CALCULATES THE CORRECTION VECTOR DX2 GIVEN 4 C G AND GX 5 C 6 I M P L I C I T REAL*8 (A-H,0-Z) 7 COHHON /MATEXP/ JEXP 8 DIMENSION G ( N C ) i GX(HC,NV), THIX (NV) 9 DIMENSION DX1(NV), DX2(NV), G X I ( 1 5 , 1 5 ) , TEHP(15,15) 10 CALL DGPROD (GX,GX,GXI,NC,NV,NC,HC,HC,15,0, 1,1) 11 CALL DINVRT (GXI,NC,15,DET,C0ND) 12 DET = DET*1.0D1**JEXP 13 WRITE (6,2) DET, COND 14 I F (DET .EQ. 0) GO TO 400 15 CALL.DGPROD (GX,GXI,TEMP,NV,NC,NC,MC,15,15, 1,0,0) 16 I F (IS .EQ. 1) GO TO 200 17 CALL DGMATV (TEMP,G,DX2,NV,NC,15) 18 I F (IS .EQ. 2) GO TO 300 19 200 CONTINUE 20 CALL DGHULT (TEMP,GX,GXI,NV,NC,NV,15,MC,15) 21 DO 260 I = 1, NV 22 DO 240 J = 1, NV 23 I F (I .EQ. J) GO TO 220 24 G X I ( I , J ) = - G X I ( I , J ) 25 GO TO 240 26 220 CONTINUE 27 G X I ( I , J ) = 1 - G X I ( I , J ) 28 240 CONTINUE 29 260 CONTINUE 30 CALL DGMATV (GXI,THIX,DX1,NV,NV,15) 31 300 RETURN 32 400 WRITE (6,1) JEXP 33 STOP 34 1 FORMAT (1X,I4,25H * * * SINGULAR MATRIX ) 35 2 FORMAT (10H DET = ,D14.7, 10H COND = , D14.7) 36 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES 60 MAIN PROGRAM 1 C COHHON /SE&RCH/ GOES HERE 2 C COHHON /SERVC/ GOES HERE 3 EQUIVALENCE (JACOB, ACOB) 4 DIMENSION JACOB (1), ISAV (200) , JSAV (200) 5 C 6 DATA NSAV /0/ 7 NCOM = (LOCF(TEMP(1))-LOCF(ACOB))/4 8 10 CONTINUE 9 IF (NSAV .EQ. 0) GO TO 60 10 DO 20 I = 1, NCOM 11 JACOB (I) = 0 12 20 CONTINUE 13 DO 40 1 = 1 , NSAV 14 J = ISAV (I) 15 JACOB (J) = JSAV(I) 16 40 CONTINUE 17 60 CONTINUE 18 CALL RSEARC 19 NSAV = 0 20 DO 80 I = 1, NCOM 21 IF (JACOB(I) .EQ.,0) GO TO 80 22 NSAV = NSAV + 1 23 ISAV (NSAV) = I 24 JSAV (NSAV) = JACOB (I) 25 80 CONTINUE 26 CALL MINMTS 27 GO TO 10 28 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES BLOCK DATA 61 1 BLOCK DATA 2 C COMMON /SEARCH/ GOES HERE 3 C COMMON /SEHVC/ GOES HERE 4 DATA TEMP /25*0.0 / 5 DATA STEMP/25*0. 9 / 6 DATA I /o / 7 DATA J /o / 8 DATA L /o / 9 DATA H /o / 10 DATA N /o / 11 . DATA NULL /1HU / 12 DATA N00 /o / 13 DATA N01 /1 / 14 DATA N02 /2 / 15 DATA N03 /3 / 16 DATA N04 / 17 DATA N05 /5 / 18 DATA N06 /6 / 19 DATA N07 /7 / 20 DATA N08 /8 / 21 DATA N09 /9 / 22 DATA N10 /10 / 23 DATA FPP5 / .5 / 24 DATA FP1 /1.0 / 25 DATA FP2 /2.0 / 26 DATA FP3 /3.0 / 27 DATA FP4 /4. 0 / 28 DATA FP5 /5.0 / 29 DATA FP6 /6.0 / 30 DATA FP7 /7.Q / 31 DATA FP8 /8.0 / 32 DATA FP9 /9.0 / 33 DATA FP10 /10.0 / 34 c DATA RPD /1712 435 7506504 512351B/ 35 c DATA DPR /17257122734064617020B/ 36 DATA XINF /100000000000.0/ 37 c 38 c 39 DATA CONEPS/ 89. 9 i 40 1 5*^ 1 / 41 DATA CONSEX/ .000001 / 42 DATA CTHA / • 5 / 43 DATA FITERR/ .000001 / 44 DATA GAHAX / 3. 0 / 45 DATA ICGM / 63 / 46 DATA ISTART/ 63 / 47 DATA IDEB / 0 / 48 DATA IDEPVR/ 25*0 / 49 DATA KASE / 1 / 50 DATA MAXITR/ 10 / 51 DATA MODEM / 1 / 52 DATA NDEPV / 0 / 53 DATA NFLAG / 0 / APPPENDIX III - CONSTRAINED TARGETING ROUTINES BLOCK DATA 62 54 DATA NINDV / 0 / 55 DATA OPT / 0.0 / 56 DATA PCTCC / .3 / 57 DATA PCTOLD/ .3 / 58 DATA PERT /25*1.0E-4 / 59 DATA PGEPS / 1.0 / 60 DATA P2HIN / 1.0 / 61 DATA SRCHM / 0.0 / 62 DATA STMINP/ .01 / 63 DATA STPHAX/ 1.0E+10 / 64 DATA D / 25*0.0 / 65 DATA ICON / 100.0 / 66 DATA HOPT / 1.0 / 67 DATA 8U / 25*1.0 / 68 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE ABT 63 1 SUBROUTINE ABT (A,B,C,L,M,N) 2 C 3 C*** PERFORMS MATRIX MULTIPLICATION 4 C C = A •* TRS (B) 5 C MATRICES ARE ASSUMED TO BE STORED IN COLUMN ORDER 6 C DIMENSION A(L,M), B (N, M) , C (L, N) 7 C 8 DIMENSION A(1), B(1), C(1) 9 C 10 DO 60 I = 1, L. 11 DO 40 J = 1, N 12 IC = I • (J-1) *L 13 C(IC) = 0.0 14 DO .20 -K = 1, fl 15 IA = (K-1) *L 16 IB = J..* (K-1)*N 17 C(IC) = C(IC) • A(IA)*B(IB) 18 20 CONTINUE 19 40 CONTINUE 20 60 CONTINUE 21 RETURN 22 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE AUXRK 64 1 SUBROUTINE AUXBK (X,F) 2 COMMON /COM/ C (12) 3 DIMENSION X (1) ,F (1) ,S (3) ,PCT (3) 4 DATA » /2.5E1,7.5E1,5.46E2/ 5 DO 20 1 = 1 , 3 6 14 = 4* (1-1) + 1 7 T = 0.0E0 8 DO 10 J = 1, 4 9 K = 14 + J 10 T = T + ABS(X(K)) 11 10 CONTINUE 12 IF (T .NE. 0) PCT(I) = ABS (X (14+1) )/T -2.0E-1 13 PCT(I) = PCT(I) *1.0E-1 14 IF (PCT(I) .GT. 0) PCT(I) = O.OEO 15 20 CONTINUE 16 DO 40 I = 2, 13 17 12 = I + 2 18 J = 12/4 19 K = 3*(I2-4*J) 20 XX = X (I) 21 IF (XX • LT. -1) XX = 1.0/ABS(XX)**3 22 IF (XX .LT. 0) XX = ABS (XX) 23 F(I) = (C(K+1) + C(K+2)*X(1) + C(K+3)*W(J) + PCT(J) ) *XX 24 C IF (F(I) .GT. 0) F (I) = 0 25 40 CONTINUE 26 RETURN 27 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE BTH 65 1. SUBROUTINE BTW (B,W,D,L,M,N) 2 C 3 C*** PERFORMS THE INDICATED MATRIX MULTIPLICATION 4 C MATRICES ARE ASSUMED TO BE STORED IN COLUMN ORDER. 5 C D = TRS (B) *B 6 C DIMENSION B(H,L), B(M,N), D (L, N) 7 C 8 DIMENSION B(1), 8(1), D(1) 9 DO 60 1=1,L 10 IM1M=M*(I-1) 11 DO 40 J=1,N 12 JM1L=L*(J-1) 13 JM1M=fl* (J-:1) 14 ID=JH1L+I 15 D(ID)=0.0 16 DO 20 K=1,M 17 IB=IM1M+K 18 IB=JM1M+K 19 D(ID)=D(ID)+B(IB)*8(IH) 20 20 CONTINUE 21 40 CONTINUE 22 60 CONTINUE 23 RETURN 24 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE BUCKET 66 1 . SUBROUTINE BUCKET(X,Y,N,XX,YY,NP) 2 C*** 3 C 4 C THIS ROUTINE REARRANGES AN ARRAY X AND THE CORRESPONDING ELEMENTS 5 C OF Y IN ASCENDING ORDER. N IS THE HUHBER OF ELEMENTS IN EACH. 6 C NP IS A POINTER TO THE FIRST ELEMENT OF Y THAT IS LESS THAN THE 7 C NEXT ELEMENT. , NP IS ZERO IF IT DOES NOT EXIST.J ORDERED VALUES 8 C ARE RETURNED IN XX AND YY. 9 C 10 C*** 11 DIMENSION X(1), Y(1), XX (1) , YY(1)jXS(10) 12 DO 10 I = 1, N .13 10 XS (I) = X(I) 14 DO 30 I = 1, N 15 NSMAL = 1 16 DO .20 J = 1, N 17 IF (XS(NSMAL) . GT. XS (J)) NSMAL = J 18 20 CONTINUE 19 XX (I) = X(NSMAL) 20 YY (I) = Y(NSMAL) 21 XS (NSMAL) = 1.E10 22 30 CONTINUE 23 NP = 0 24 DO 40 I = 2, N 25 IF (YY(I-1) .LE. YY (I) ) GO TO 50 26 40 CONTINUE 27 RETURN 28 50 NP = 1-1 29 RETURN 30 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE CGM 67 1 SUBROUTINE CGM (F0,G,GHAG) 2 C THIS SUBROUTINE CALCULATES A DU BI THE CONJUGATE GRADIENT METHOD 3 C WITH AN INITIAL STEEPEST DESCENT STEP WHEN ICGH IS NONZERO 4 C 5 C COMMON /SEARCH/ GOES HERE 6 C COMMON /SERVC/ GOES HERE 7 EQUIVALENCE (GSQ,TEHP(1) ) , (CON,TEMP (2) ) 8 C 9 DIMENSION G(1) 10 GSQ=GMAG**2 11 CON=GSQ/GSQOLD 12 IF (ICGM.EQ.0) GO TO 40 13 C 14 C INITIALIZATION PASS 15 C 16 L=0 17 DO 20 I=1,NINDV 18 S(I)=-G(I) 19 DU(I)=S(I) 20 IF (IDAV.EQ.O) GO TO 20 21 L=L+I 22 H(L) = FP1 23 GP(I) = G(I) 24 20 CONTINUE 25 GO TO 100 26 C 27 C CONJUGATE STEPS 28 C 29 40 CONTINUE 30 IF (IDAV.EQ.O) GO TO 60 31 CALL DGM (G) 32 GO TO 100 33 60 CONTINUE 34 DO 80 I=1,NINDV 35 S(I)=CON*S(I)-G(I) 36 DD(I)=S(I) 37 80 CONTINUE 38 100 CONTINUE 39 GSQOLD=GSQ 40 ICGM=0 41 RETURN 42 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES 68 SUBROUTINE COMBIN 1 SUBROUTINE COMBIN 2 C COMMON /SEARCH/ GOES HERE 3 IF (IND(1) .NE. 0) GO TO 40 4 M = IND(2) 5 N = IND(3) 6 DO 20 I = 1 , N 7 IND(I) = I 8 20 CONTINUE 9 1 = 1 10 K = N 11 GO TO 120 12 40 CONTINUE 13 IF (K .EQ. 0) GO TO 100 14 L = M + K - N 15 IF (IND(K) .LT. L) GO TO 60 16 K = K - 1 17 GO TO 40 18 60 CONTINUE 19 IND(K) •= IND(K) +1 20 1 = 1 + 1 21 80 CONTINUE 22 IF (K .EQ. N) GO TO 120 23 IND(K+1) = IND(K) • 1 24 K = K + 1 25 GO TO 80 26 100 CONTINUE 27 IND(1) = -1 28 120 CONTINUE 29 RETURN 30 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE CUBMIN 69 1 SUBROUTINE CUBMIN (A, XMIN,YMIN) 2 C 3 C INPUTS 4 C A=ARRAY OF CUBIC POLYNOMIAL COEFFICIENTS. A(1) IS CONST 5 C TERM, A (2) IS LINEAR TERM, A (3) IS QUADRATIC TERM, AND 6 C A(4) IS CUBIC TERM. 7 C 8 C OUTPUTS 9 C XMIN=ABSCISSA VALUE OF MINIMUM OF CUBIC 10 C YMIN=ORDINATE VALUE OF MINIMUM OF CUBIC 11 C 12 1 FORMAT( /,47X,36HCUBIC DEGENERATES TO A STRAIGHT LINE,/) 13 2 FORMAT( /,55X,20HCUBIC HAS NO EXTREMA,/) 14 3 FORMAT( /,41X,41HCUBIC DEGENERATES TO A QUADRATIC WITH NO 15 1 , 7HMINIMUM,/) 16 C 17 DIMENSION 18 1 A(4) 19 C 20 DATA 21 1 R2 /2. / 22 2 ,R3 /3. / 23 3 ,TP10 /1.0E+10 / 24 C 25 IF (A (4)) 10,20,10 26 20 CONTINUE 27 IF (A (3)) 30,40,50 28 40 CONTINUE 29 WRITE (6,1) 30 GO TO 60 31 30 CONTINUE 32 WRITE (6, 3) 33 GO TO 60 34 50 CONTINUE 35 XMIN=-ft(2)/(R2*A(3)) 36 GO TO 70 37 10 CONTINUE 38 THRA4=R3*A(4) 39 DISC=A (3) **2-A (2) *THRA4 40 IF(DISC) 80,80,90 41. 80 CONTINUE 42 C WRITE (6,2) 43 60 CONTINUE 44 XMIN=TP10 45 YMIN=TP10 46 GO TO 100 47 90 CONTINUE 48 DISC=SQRT(DISC) 49 XMIN=(-A(3)+DISC)/THRA4 APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE CUBHIN 70 50 70 CONTINUE 51 YHIN=A (1) +XMIN* (A (2) +XMIN* (A (3) +XHIN*A (4) ) ) 52 100 CONTINUE 53 RETURN 54 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE DELTU 71 1 SUBROUTINE DELTU 2 C 3 C DELTU CALCULATES THE DIRECTION OF SEARCH BY THE METHOD SPECIFIED 4 C BY SRCHM, SCALES THE VECTOR AND CALLS ULIMIT AND TEST 5 C 6 C COMMON /SEARCH/ GOES HERE 7 INTEGER SRCHM 8 C 9 CALL BUCAL 10 CALL GMAG 11 C 12 C BRANCH ON SRCHM 13 C 14 GO TO (20,40,60,80)SRCHM 15 C CONJUGATE GRADIENT METHOD 16 20 CONTINUE 17 ICGM=63 18 CALL CGM (P2N0M,G2,G2MAG) 19 GO TO 100 20 C STEEPEST DESCENT METHOD 21 40 CONTINUE 22 CALL SDM 23 GO TO 100 24 C CGM OF G1 25 60 CONTINUE 26 ICGM=63 27 CALL CGM (P1N0M,G1,G1MAG) 28 GO TO 100 29 C PROJECTED GRADIENT METHOD 30 80 CALL PGM 31 100 CONTINUE 32 CALL UNITDU 33 IF (SRCHM .EQ. 4) CALL GABDD 34 RETURN 35 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE DGM 72 1 SUBROUTINE DGM (G) 2 C 3 C*** COMPUTES DO BY DAVIDONS DEFLECTED GRADIENT METHOD 4 C 5 C 6 C COMMON /SEARCH/ GOES HERE 7 C COMMON /SERVC/ GOES HERE 8 C 9 EQUIVALENCE (P,DU) , (PP,S) 10 DIMENSION P(25), PP (25) , G(25) 11 DATA ZERO,ONE/0. 0,-1.0/ 12 N = NINDV 13 M=N*(N+1)./2 14 C COMPOTE NEB DEFLECTION MATRIX 15 DO 20 J=1,N 16 DG(J)=G(J)-GP(J) 17 20 CONTINUE 18 DO 100 J=1,B 19 HDG(J)=ZERO 20 DO 100 L=1,N 21 IF (J-L) 60,60,40 22 40 JLA=L+J*(J-1)/2 23 GO TO 80 24 60 JLA=J+L* (L-1)/2 25 80 HDG(J)=HDG(J) + H(JLA) *DG (L) 26 100 CONTINUE 27 PPTDG=ZERO 28 DGTHDG-ZERO 29 DO 120 J=1,N 30 PPTDG=PPTDG+PP (J) *DG (J) 31 DGTHDG=DGTHDG+DG(J) *HDG (J) 32 120 CONTINUE 33 TEMP 1=GAMASS/PPTDG 34 TEMP2=ONE/DGTHDG 35 DO 140 L=1,N 36 DO 140 J=1,L 37 JLA=J+L*(L-1)/2 38 A=PP(J) *PP(L) *TEMP1 39 B=HDG(J) *HDG(L) *TEMP2 40 H (JLA) =H (JLA) +A-B 41 140 CONTINUE 42 C COMPUTE NEW SEARCH DIRECTION 43 DO 220 J=1,N 44 P(J)=ZERO 45 DO 220 L=1,N 46 IF • (J-L) 180,180,160 47 160 JLA=L+J* (J-1)/2 48 GO TO 200 49 180 JLA=J+L*(L-1)/2 50 200 P(J)=P(J)-H(JLA)*G(L) 51 220 CONTINUE 52 DO 240 J=1,N APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE DGH 73 53 PP(J)=P(J) 54 GP(J)=G(J) 55 240 CONTINUE 56 RETURN 57 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE DIGDIF 74 1 SUBROUTINE DIGDIF(M,N,NDIF) 2 DIMENSION MTAB (8) 3 DATA MTAB /9999999,999999,99999,9999,999,99, 9,0/ 4 L = LXOR(M,N) 5 DO 10 I = 1, 8 6 IF (L .GT. MTAB (I)) GO TO 20 7 10 CONTINUE 8 1 = 9 9 20 NDIF = 9 - 1 10 WRITE (6,1) H,N,L,I 11 1 FORMAT (10H H,N,L,I , 4 (5X,Z8) ) 12 RETURN 13 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE DISPLY 75 1 SUBROUTINE DISPLY 2 CC COMMON /SEARCH/ GOES HERE 3 CC COMMON /SERVC/ GOES HERE 4 C EQUIVALENCE (IR1,I) 5 CC 6 C DIMENSION MSG (3) 7 C DATA IPP/1523070000Q10000OO00B/ 8 C DATA LRAP1/0/ 9 C IF (LRAP1.NE.0) GO TO 20 10 C LRAP1=2-L0CF(MSG) 11 C IR1=LOCF(MSG) 12 C IPP=IPP+IR1 13 C20 CONTINUE 14 C TEMP(1)=ABS(P1N0M) 15 C ENCODE (12,80,STEHP(1) ) TEMP{1) 16 C ENCODE (12,80,TEHP(1)) P2NOM 17 C IF (P2NOM.EQ.Q.0) P2NOH=1.0E-10 18 C IR1 = 0 19 C IF (P2NOM . NE. 0.0) IR1 = ALOG10(ABS(P2NOH)) 20 C 1=1H+ 21 C IF (P2NOM.GE. 1.0) GO TO 40 22 C I=1H-23 C IR1=1-IR1 24 C40 CONTINUE 25 C ENCODE (28#100rMSG) NSTEP,STEMP(1)#TEMP(1), I,IR1#CTHA 26 C60 CONTINUE 27 C IF (MSG (LRAP1) . NE. 0) GO TO 60 28 C MSG (LRAP1) = IPP 29 RETURN 30 C . 31 C80 FORMAT (E12.6) 32 C100 FORMAT (12,1X/A8,1X,A4,1HE,A1#I1,1X,F8.5) 33 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE DROP 76 1 SUBROUTINE DROP 2 C COMMON /SEARCH/ GOES HERE 3 C COMMON /SERVC/ GOES HERE 4 C 5 IF (NAC .GT.NINDV) GO TO 40 6 20 CONTINUE 7 CALL REVISE 8 IF (IDEB . NE. 0) WRITE (6,3) ICD 9 IF (ICD .EQ.,0) RETURN 10 GO TO 20 11 C 12 C MORE ACTIVE CONSTRAINTS THAN CONTROLS 13 C 14 40 CONTINUE 15 DP1 = -1.0E+10 16 IND(1) = 0 17 IND(2) = NAC - NEQC 18 IND(3) = NINDV - NEQC 19 IF (IND(2) *IND(3) .LE. 0) GO TO 300 20 DO 60 I = 1, NDEPV 21 ISTC(I) = ITC(I) 22 60 CONTINUE 23 C 24 C MAKE ONE ITERATION PER COMBINATION 25 C 26 80 CONTINUE 27 CALL COMBIN 28 IF (IND(1) .LT. 0) GO TO 200 29 J = 0 30 K = 1 31 DO 100 I = 1, NDEPV 32 ITC(I) =0 33 IF (ISTC(I) .EQ. 0) GO TO 100 34 J = J + 1 35 IF (J .NE. IND(K)) GO TO 100 36 K = K + 1 37 ITC(I) • = 63 38 100 CONTINUE 39 NTC = NINDV-NEQC 40 WRITE (6,5) (ITC (I) ,1=1,NDEPV) 41 CALL UPDATS 42 CALL REVISE 43 IF (ICD .EQ. 0) GO TO 400 44 WRITE (6,7) (ITC (I) ,1=1,NDEPV) 45 DO 140 J = 1, NDEPV 46 IF (ISTC(J) .EQ- ,0) GO TO 140 47 IF (ITC(J) . NE. 0) GO TO 140 48 DO 120 I = 1, NINDV 49 K = J + (1-1) *NDEPV 50 TEMP (I) = ACOB (K) 51 120 CONTINUE 52 CALL MAT PI (PG1,TEMP (1) ,STEMP (1) ,N01, NINDV , N01) APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE DROP 77 53 WRITE (6,6) STEHP(1) 54 IF (STEHP (1) ) 140,140,80 55 140 CONTINUE 56 TEHP(1)= SP (PG1,G1, NINDV)/PG1 MAG 57 IF (IDEB .NE. ,0) WRITE (6,4) TEMP (1) , (ITC (I) ,1=1, NDEPV) 58 IF (TEMP(1) .LT. DP1) GO TO 80 59 DP1 = TEMP(1) 60 DO 150 1 = 1 , NDEPV 61 IBTC(I) = ITC (I) 62 150 CONTINUE 63 GO TO 80 64 C 65 C FEASIBLE DIRECTION FOUND 66 C 67 200 CONTINUE 68 DO 220 I = 1, NDEPV 69 ITC (I) = IBTC(I) 70 220 CONTINUE 71 CALL UPDATS 72 RETURN 73 C 74 C ERROR MESSAGES 75 C 76 300 CONTINUE 77 WRITE (6,1) 78 GO TO 500 79 400 CONTINUE 80 WRITE (6,2) 81 500 CONTINUE 82 NFLAG = -1 83 RETURN 84 1 . FORMAT (44H-*** MORE EQUALITY CONSTRAINTS THAN CONTROLS) 85 2 FORMAT (49H-*** PROBLEM SOLVED === NAC IS GREATER THAN NINDV) 86 3 FORMAT (10H-CONTRAINT,I2,8H DROPPED) 87 4 FORMAT (7H-DP1 = ,E21.14,7H ITC = ,25(11, 1X)) 88 5 FORMAT (20H-COMBINATION OF ITC ,25 (11,1X)) 89 6 FORMAT (37H-NEGATIVE DU DOTTED WITH CONSTRAIN = ,E21.14) 90 7 FORMAT (20H-REVISION OF ITC ,25(I1,1X)) 91 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE ERROR 78 1 SUBROUTINE ERROR (I,J) 2 C***ERR0R 3 C* ERROR-PRINT ERR/R MESSAGE AND DETERMIN IF FATAL 4 C* 5 C DATA MSK1 /77000000000000000000B/ 6 C DATA IBK1 /55000000000000000000B/ 7 CC . 8 10 CONTINUE 9 CALL PAGER (3) 10 WRITE (6,100) I,J 11 20 CONTINUE 12 STOP 13 RETURN 14 100 FORMAT (10X,10H**ERR0R** ,2X,A10,A10 ) 15 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE FGAMA 79 1 SUBROUTINE EGAMA (IS) 2 C 3 C*** DETERMINES THE P1 AND P2 ASSOCIATED SITH A PARTICULAR GAMA IN THE 4 C. DIRECTION OF SEARCH BY CHANGING THE CONTROLS . U = DO • GAMHA*DU 5 C 6 C COMMON /SEARCH/ GOES HERE 7 C COMMON /SERVC/ GOES HERE 8 C 9 GAMAS=STEP(IS) 10 CALL TRAJ 11 . P1TRY (IS)=P1 12 P2TRY(IS)=P2 13 IF (INTRI1 .EQ. 0) RETURN 14 C 15 C ESTIMATE NET PERFORMANCE 16 C 17 IF (NAC .EQ. 0) GO TO 15 18 CALL HATPY (EHAP,EA,TEHP (1) ,NINDV,NAC, 1), 19 P1TRY (IS) =P1TRY (IS) +SP (G1,TEMP (1) ^ HINDV) 20 15 CONTINUE 21 JS=IS-1 22 SAVIT(1,JS)=P1 23 SAVIT(2#JS)=P2 24 DO 20 1=1,NDEPV 25 SAVIT(I+2,JS)=E(I) 26 20 CONTINUE 27 RETURN 28 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES 80 SUBROUTINE FOPMIN 0 1 SUBROUTINE FOPMIN (X,.Y, XHIFJ,YMIH, IERR) 2 C 3 C INPUTS 4 C X=ARRAY OF ABSCISSA VALUES OF FOUR SAMPLE POINTS 5 C . Y=ARRAY OF ORDINATE VALUES OF FOUR SAMPLE POINTS 6 C 7 C OUTPUTS 8 C XMIN=ABSCISSA VALUE OF MINIMUM OF CUBIC FITTED THROUGH 9 C XMIN=FOUR SAMPLE POINTS 10 C YMIN=ORDINATE VALUE OF MINIMUM OF ABOVE CUBIC 11. C IERR=FLAG BHOSE NON-ZERO VALUE INDICATES THAT TWO OF THE 12 C . GIVEN X VALUES ARE IDENTICAL 13 C 14 DIMENSION 15 1 X(4) ,Y<4) ,A(4) 16 C 17 DATA 18 1 NO /0 / 19 2 ,N1 /1 / 20 C 21 IERR=N0 22 B21=X(2)-X(1) 23 B31=X(3)-X(1) 24 B32=X(3)-X(2) 25 B41=X(4)-X(1) 26 B42=X(4)-X(2) 27 B43=X (4)-X (3) 28 TEMP1=B21*B31*B41 29 IF (TEMPI) 10,20,10 30 10 CONTINUE 31 F1=Y (1)/TEMPI 32 TEMP1=B21*B32*B42 33 IF (TEMPI) 30,20,30 34 30 CONTINUE 35 F2=-Y(2)/TEMP1 36 TEMP1=B31*B32*B43 37 IF (TEMPI) 40,20,40 38 40 CONTINUE 39 F3=Y (3)/TEMPI 40 TEHP1=B41*B42*B43 41 F4=-Y(4)/TEMP1 42 GO TO 50 43 20 CONTINUE 44 IERR=N1 45 GO TO 60 46 50 CONTINUE 47 A(4)=-(F1+F2+F3+F4) 48 A (3) = (X (2) +X (3) +X (4) ) *F1 + (X (1) +X (3) + APPPEHDIX III - CONSTRAINED TARGETING ROUTINES 81 SOBBOOTINE FOPMIN X(4))*F2 +(X(1)+X(2)+X(4))*F 49 1 • (X (1) +X (2) + X(3))*F4 50 B21=X(2) *X(1) 51 B31=X(3) *X(1) 52 B32=X(3) *X{2) 53 B41=X(4) *X(1) 54 B42=X(4) *X(2) 55 B43=X(4)*X(3) 56 A(2)=-({B32+B42+B43)*F1. +(B31 + B41+B43)*F2 • (B21 + B41 + B42) *F3 57 1 + (B21 + B31 + B32) *F4) 58 A (1) =Y (1) - (X (1) * (A (2) +X (1) * (A (3) + X(1)*A(4)))) 59 CALL CUBHIN (A, XMIN, JTMIN) 60 60 CONTINUE 61 RETURN 62 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES 82 SUBROUTINE GABDD 1 SUBROUTINE GABDD 2 C COMMON /SEARCH/ GOES HERE 3 C . COMMON /SERVC/ GOES HERE 4 C 5 STPMAX = 1.0E+10 6 DO 40 I = 1, NDEPV 7 IF (IDEPVR (I) .EQ. 0) GO TO 40 8 IF (ENOM(I) .LE. 1.0) GO TO 40 9 DO 20 J = 1, NINDV 10 K = I + (J-1)*NDEPV 11 TEMP (J) = ACOB (K) 12 20 CONTINUE 13 CALL HATPY (TEMP (1) ,DU,STEMP (1) , 1, NINDV, 1) 14 IF (STEMP (1) . GE. 0) GO TO 40 15 TEMP (1).= -ENOM (I)/STEMP (1) 16 IF (TEMP(1) .GT. STPMAX) GO TO 40 17 STPMAX = TEMP(1) 18 40 CONTINUE 19 RETURN 20 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE GENHIN 83 1 SUBROUTINE GENHIN (X,Y,DYDX1,EPSA,YES,HIN) 2 C 3 C*** FINDS A MINIMUM TO THE RIGHT OF X(1), BUT TO THE LEFT EPSA(2). , 4 C TRIAL STEPS ARE STORED IN ORDERED PAIRS (X, Y) . DYDX1 IS THE SLOPE 5 C OF THE FUNCTION AT X(1)., ZERO SLOPE IS IGNORED IN THE ESTIMATION 6 C PROCESS. USER MAY PRESET ANY ELEMENTS OF THE X AND Y ARRAYS., 7 C YES CONTAINS THE PREDICTED VALUES OF THE FUNCTION BASED UPON THE 8 C CURVE FIT. IF ANY PREDICTED VALUE IS LESS THAN EPSA (4) PERCENTAGE 9 C FROM THE ACCUAL VALUE, THE HINIMUZATION PROCESS IS TERMINATED., 10 C IF ANY TWO CONSECUTIVE X VALUES ARE LESS THAN EPSA(3) PERCENTAGE 11 C APART, CONVERGENCE IS ASSUMED. MIN IS SET TO THE NUMBER OF THE 12 C ORDERED PAIR THAT COMES CLOSEST TO THE MINIMUM. A NEGATIVE MIN 13 C INDICATES THAT A WELL DEFINED MINIMUM DID NOT EXIST. A ZERO MIN 14 C INDICATES THAT NO MINIMUM EXISTS. , EPSA (2) IS THE MAXIMUM PERCENT 15 C THAT X IS ALLOWED TO BE REDUCED BY IF THE FUNCTION INCREASES 16 C . DRASTICALLY. 17 C 18 DIMENSION X (6) , Y (6) , XX(5), YY(5), EPSA(4), ^ YES (4) 19 DO 20 1=1,4 20 YES (I) =0.0 21 20 CONTINUE 22 IF (DYDX1) 40,360,440 23 CSLOPE NEGATIVE 24 CFILL IN NEEDED VALUES FOR FIRST FIT 25 40 IF (Y(1) -EQ.0.0) CALL .FGAMA (1) 26 IF (Y(2) .NE.0.0) GO TO 60 27 IF (X(2) .EQ.0.0) X(2)=1.0 28 CALL FGAMA (2) 29 GO TO 80 30 60 IF <Y{3) .NE.O.Q) GO TO 120 31 CTWO POINT ONE SLOPE 32 80 CALL TPOSM (X,Y,DYDX1,X (3) ,YES) 33 IF (X(3) .LT.EPSA(1)*X(2)) X (3) =EPSA (1) *X (2) 34 IF (X(3) .GT.EPSA (2) ) X(3)=EPSA(2) 35 IF (ABS((X(3)-X(2))/X(2)).GT.EPSA(3)) GOTO 10P 36 MAX=2 37 GO TO 300 38 100 CALL FGAMA (3) APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE GENMIN 84 39 CCHECK QUADRATIC FIT 40 IF (Y{3) .EQ.0.0) Y(3)=1.0E-10 41 IF (ABS ( (YES (1).—¥ (3) )/Y (3) ).GE.EPSA (4) ) GO TO 120 42 HIN=3 43 GO TO 460 44 CTHREE POINT ONE SLOPE 45 120 CALL THPOSH (X, Y,D¥DX1#X (4) #YES (2) ,ISER) 46 140 CONTINUE 47 IF (X(4) .LT.EPSA (1) *X(3) ) X (4) =EPSA (1) *X (3) 48 IF (X(4) .GT.EPSA (2) ) X(4)=EPSA(2) 49 IF (ABS ( (X (4) -X (3) ) /X (3) ) . LE. EPSA (3) ) GO TO 150 50 IF (ABS ( (X (4) -X (.2) ) /X (2)) .GT. EPSA (3) ) GO TO 160 51 150 MAX=3 52 GO TO 300 53 160 CALL FGAMA (4) 54 CCHECK CUBIC FIT 55 IF (Y (4) .EQ.0.0) Y(4) = 1.0E-10 56 IF (ABS((YES(2)~Y(4))/Y(4)). GE.EPSA (4)) GO TO 180 57 MIN=4 58 GO TO 460 59 CFIND THREE BEST POINTS FOR NEXT FIT 60 180 CALL BUCKET (X,¥#4#XX#YY,NP) 61 IF (NP-GT. 1) GO TO 220 62 200 MAX=4 63 GO TO 300 64 220 NP=NP-1 65 CTHREE POINT 66 CALL THPH (XX (NP) , YY (NP) ,X (5) , YES (3) ) 67 IF (X (5) .XT. EPSA (1) *X (4) ) X (5) =EPSA (1) *X (4) 68 IF (X(5) .GT.EPSA (2) ) X(5)=EPSA(2) 69 IF (ABS ((X (5) -X (4) ) /X (4) ) . LE. EPSA (3) ) GOTO 200 70 IF (ABS((X(5)-X(3))/X(3)) .LE.EPSA(3)) GOTO 200 71 IF (ABS((X(5)-X(2))/X(2)).LE.EPSA(3)) GOTO 200 72 CALL FGAMA (5) 73 CCHECK QUADRATIC (NO SLOPE) FIT 74 IF (Y (5) .EQ.0.0) Y(5) = 1.0E-10 75 IF (ABS((¥ES(3)-Y(5))/Y(5)) .GE.EPSA (4)) GO TO 240 76 MIN=5 77 GO TO 460 78 240 XX (NP+3)=X(5) 79 YY (NP+3) =Y (5) 80 CFOUR POINT CUBIC FIT 81. CALL FOPMIN (XX (NP) ,YY (NP) , X (6) , YES (4) ,ISER) 82 IF (X(6) .LT.EPSA (1) *X(5)) X (6) =EPSA (1) *X (5) 83 IF (X (6) .GT. EPSA (2) ) X (6) =EPSA (2) 84 IF (ABS((X(6)-X(5) )/X(5) ) .GT.EPSA(3) ) GOTO APPPENDIX III — CONSTRAINED TARGETING ROUTINES 85 SUBROUTINE GENHIN 260 85 MAX=5 86 GO TO 300 87 260 CONTINUE 88 CALL FGAMA (6) 89 C CHECK CUBIC (NO SLOPE) FIT 90 IF (Y(6)-EQ.0.0) >Y (6) =1. OE-10 91 IF (ABS((YES(4)-Y(6))/Y(6)) .GE.EPSA(4)) GO TO 280 92 MIN=6 93 GO TO 460 94 280 MAX=6 95 300 MIN=1 96 CFIND BEST POINT AND RETURN 97 320 HINSV=MIN 98 DO 340 1=1,MAX 99 IF (Y(MIN) .GT.Y(I)) MIH=I 100 340 CONTINUE 101 IF (MIN.EQ.1) GO TO 440 102 IF (MINSV.NE.MIN) MIN=-MIN 103 RETURN 104 CSLOPE OF ZERO INDCATES ONLY POINTS ARE TO BE USED FOR THE CURVE FIT 105 360 IF (1(1) .NE.0.0) GO TO 380 106 CALL FGAMA (1) 107 380 IF (Y(2) .NE.0.0) GO TO 400 108 IF (X(2) .EQ.0.0) X(2)=.?5 109 CALL FGAMA (2) 110 400 IF (Y (3) .NE.0.0) GO TO 420 111 IF (X(3) .EQ.0-0) X (3) =1.25 112 CALL FGAMA (3) 113 420 CALL THPM (X, Y,X (4) ,IES (2) ) 114 GO TO 140 115 CSLOPE POSITIVE 116 440 MIN=0 117 RETURN 118 460 MAX=MIN 119 GO TO 320 120 END APPPENDIX III - CONSTRAINED TARGETING ROOTINES SUBROUTINE GMAG 86 1 SUBROUTINE GHAG 2 C 3 C*** CALCULATES THE MAGNITUDE OF THE GRADIENTS 4 C 5 C COMMON /SEARCH/ GOES HERE 6 C COMMON /SERVC/ GOES HERE 7 C 8 INTEGER SRCHH 9 IF (OPT.EQ.O.OH.SRCHM.EQ.1) GO TO 20 10 G1MAG=SQRT(SP(G1,G1 ,NINDV) ) 11 IF (G1MAG.EQ.O.0) CALL ERROR (10H G1MAG=0, , 10HCK INPUTS ) 12 C 13 20 CONTINUE 14 IF (NDEPV.EQ.O) GO TO 40 15 \ G2MAG=SQRT(SP(G2,G2,NINDV) ) 16 40 CONTINUE 17 RETURN 18 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES 87 SUBROUTINE GHAD 1 SUBROUTINE GRAD 2 C 3 C GRAD CALCULATES THE GRADIENTS TO EACH OF THE TARGETS AND TO THE 4 C OPTIMIZATION INDEX WITH RESPECT TO THE CONTROLS 5 C 6 C COMMON /SEARCH/ GOES HERE 7 C COMMON /SERVC/ GOES HERE 8 INTEGER SRCHM 9 EQUIVALENCE (ISUB,I) 10 C 11 GAMAS=1.0 12 DO 20 I=1#NINDV 13 DU (I) =0.0 14 20 CONTINUE 15 DO 140 K=1#NINDV 16 DU (K)=PERT(K) 17 IF (NDEPV.EQ.O) GO TO 60 18 DO 40 L=1,NDEPV 19 E(L)=ENOH(L) 20 40 CONTINUE 21 60 CONTINUE 22 P1=P1N0H 23 C . 24 KS = K 25 CALL TRAJ 26 K = KS 27 C 28 DU (K)=0.0 29 G2 (K)=0.0 30 IF (OPT.EQ.0.OR.SRCHM.EQ.1) GO TO 80 31 G1 (K) = (P1-P1N0H)/PERT(K) 32 CALL PAD (P1NOM#P1#0) 33 80 CONTINUE 34 IF (NDEPV.EQ.O) GO TO 120 35 DO 100 L=1,NDEPV 36 ISUB=L+(K-1)*NDEPV 37 CALL PAD (ENOH (L) ,E (L) , 0) 38 ACOB(ISOB) = (E (L) -ENOM (L) ) /PERT (K) 39 G2(K) = G2(K) • FP2*AC0B (ISUB) *ENOM (L) 40 100 CONTINUE 41 . 120 CONTINUE 42 CALL PAD (PERT (K) ,U (K) , N01) 43 140 CONTINUE 44 RETURN 45 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE INVH 88 1 SUBROUTINE INVH (A,N) 2 C 3 C INPUT 4 C A (N, N) MATRIX 5 C N = DIMENSIONS OF A 6 C OUTPUT 7 C, A (N,N) MATRIX ***** A (OUTPUT.) • = INVERSE OF A(INPUT) 8 C 9 DIMENSION A(1) ,XB(25) ,RTEST(25) ,IX(25) , IY(25) 10 CALL ZERO (RTEST,25,1) 11 CALL ZERO (XB,25,1) 12 C 13 C MAIN LOOP 14 C 15 DO 100 1=1,N 16 S1 = 0. 17 DO 40 JJ=1,N 18 IF (XB(JJ)) 40,1,40 19 1 LA = (JJ-1) *N 20 DO .4 J=1,N 21 IF (RTEST(J)) 4,2,4 22 2 NB = LA+J 23 S2 = ABS (A (NB) ) 24 IF (S1-S2) 3,4,4 25 3 S1 = S2 26 LC = JJ 27 NC = J 28 NG = NB 29 4 CONTINUE 30 40 CONTINUE 31 IF (S1) 5,5,6 32 5 CALL ZERO (A,N,N) 33 WRITE (6,200) 34 200 FOBMAT (45HOINVERSION OF A SINGULAR MATRIX WAS ATTEMPTED) 35 RETURN 36 6 IX (I) = NC 37 IY (I) = LC 38 RTEST(NC) = 10. 39 XB(LC) = A (NG) 40 C 41 C CLOBBER LOOP 42 C 43 NAA = (LC-1) *N 44 DO 8 K=1,N 45 IF (NC-K) 9,8,9 46 9 DO 13 KL=1,N 47 IF (KL-LC) 131,13,131 48 131 NI = (KL-1)*N 49 NJ = NAA+K 50 NF = NI+K APPPENDIX III — CONSTRAINED TARGETING ROUTINES 89 SUBROUTINE INVM 51 NH = NI+NC 52 A (NF) = A(NF)-A(NJ)*A(NH)/XB(LC) 53 13 CONTINUE 54 8 CONTINUE 55 DO 11 .KK=1,N 56 NE = NAA+KK 57 11 A(NE) = -A(NE)/XB(LC) * 58 A(NG) = 1. 59 C 60 C END OF CLOBBER LOOP 61 C 62 100 CONTINUE 63 C 64 C END OF MAIN LOOPS, START OF POSITIONING LOOPS 65 C 66 C . (500 LOOP IS COLUMN REPOSITIONING) 67 C 68 DO 500 1=1,N 69 DO 501 J=1,N 70 NT = (IY(J)-1) *N+I 71 501 RTEST (J) = A (NY) 72 DO 500 J=1,N 73 NY = (IX(J)-1) *N+I 74 500 A(NY) = RTEST(J) 75 C 76 C (600 LOOP IS ROD REPOSITIONING) 77 C 78 DO 600 1=1,N 79 NX = (1-1) *N 80 DO 601 J=1,N 81 NY = NX+IX(J) 82 NZ = IY(J) 83 601 RTEST (NZ) = A (NY) 84 DO 600 J=1,N 85 NZ = NX+J 86 600 A(NZ) = RTEST (J) 87 C 88 C END OF POSITIONING LOOPS ***** START COMPUTING INVERSE 89 C 90 DO 800 1=1,N 91 D = XB (I) 92 DO 800 J=1,N 93 NZ = (J-1) *N+I 94 800 -A (HZ). » A(NZ)/D 95 C 96 RETURN 97 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE ITERO 90 1 SUBROUTINE ITERO 2 C 3 C PROGRAM ITERO IS THE MAI 4 C THE RESULTS OF THE ITERATION 5 C 6 C COMMON /SEARCH/ GOES HERE 7 C COMMON /SERVC/ GOES HERE 8 INTEGER SRCHM 9 C 10 C LABLES DATED IN 11 C 12 DATA LCTHA /4H CTH/ 13 DATA LCP /4H CP// 14 DATA LDP1DS/4H DP1/ . 15 DATA LDP2DS/4H DP2/ 16 DATA LDU /4H DU / 17 DATA LDUMAG/4H DUM/ 18 DATA LB /4H E(1/ 19 DATA LGAMA /4H GAM/ 20 DATA LGAMAS/4H GAM/ , 21 DATA LG1 /4H G1(/ 22 DATA LG1MAG/4H G1M/ 23 DATA LG2 /4H G2(/ 24 DATA LG2MAG/4H G2H/ 25 DATA LIAC /4H IAC/ 26 DATA LNAC /4H NAC/ 27 DATA LPCTCC/4H PCT/ 28 DATA LPERT /4H PER/ 29 DATA LPG1 /4H PG1/ 30 DATA LPG1M /4H PG1/ 31 DATA LP1 /4H P1 / . 32 DATA LP1TRY/4H P1T/ 33 DATA LP2 /4H P2 / 34 DATA LP2TRY/4H P2T/ 35 DATA LRATIO/4H RAT/ 36 DATA LSMAT /4H SMA/ 37 DATA LSTPMX/4H STP/ 38 DATA LU /4H U (1/ 39 DATA LUMAG /4H UMA/ 40 DATA L«P1 /4H 8P1/ 41 DATA LBVEC /4H 8VE/ 42 DATA L8U /4H iU{/ 43 DATA LYPRED/4H YPR/ 44 C 45 C OUTPUT BLOCK 46 C 47 IF (NOMF.NE.O) GO TO 20 48 L=(NINDV+5)/6 49 L=3*L 50 M= (NDEP?+5)/6 51 M=2*M 52 I = L+M+11 53 CALL PAGER(I) APPPENDIX III — CONSTRAINED TARGETING ROUTINES 91 SUBROUTINE ITERO 54 WRITE (6,620) 55 GO TO 300 56 20 CONTINUE 57 CALL PAGER (1000) 58 IF (NSTEP.NE.1) GO TO 40 59 WRITE (6,480) 60 GO TO 300 61 40 CONTINUE 62 M=NSTEP-1 63 WRITE (6,500) M 64 80 CONTINUE 65 TEMP (1)=SECONS 66 CALL TIME (1,0,IJKLMN) 67 SECONS = IJKLMN 68 TEMP (1) =SECONS—TEMP (1) 69 WRITE (6,560) LCP,TEHP(1) 70 IF (IDEB.EQ.O) GO TO 100 71 WRITE (6,560) XPERT, (PERT (I) ,1=1,NINDV) 72 100 CONTINUE 73 IF (NDEPV.EQ.O) GO TO 140 74 DO 120 J=1,NDEPV 75 WRITE (6,560) LSMAT, (ACOB ( J+I*NDEPV-NDEPV), I=1,NINDV) 76 120 CONTINUE 77 140 CONTINUE 78 IF (OPT.EQ.0.0.OR.SRCHM.EQ.1) GO TO 160 79 WRITE (6,560) LG1, (G1 (I) ,1=1 ,NINDV) 80 WRITE (6,560) LG1HAG,G1MAG 81 160 IF (NDEPV.EQ.O) GO TO 180 82 WRITE (6,560) LG2, (G2 (I),1=1, NINDV) 83 WRITE (6,560) LG2MAG,G2MAG 84 180 IF (OPT.EQ.0.0.OR.SRCHH.EQ.1) GO TO 200 85 WRITE (6,560) LPG1, (PG1(I) ,1=1,NINDV) 86 WRITE (6,560) LPG1H,PG1MAG 87 200 CONTINUE 88 IF (IDEB.EQ.1.AND.NSTEP.EQ.2) WRITE (6,560) LWVEC, (WU (I),1=1,NINDV 89 1) 90 IF (ITERF .EQ..0) GO TO 280 91 IF (NEQC .EQ. NDEPV) GO TO 220 92 J = JAC(1) 93 WRITE (6,580) LNAC, J 94 WRITE (6,580) XIAC, (JAC (1*1) ,I=1,J) 95 220 CONTINUE 96 IF (IDEB .EQ. 1) WRITE (6,560) LRATIO, (RATIO (I) ,1=1, 4) 97 TEMP (1)=ARCOS (CTHA) *DPR 98 WRITE (6,560) LCTHA,TEMP(1) 99 WRITE (6,560) LDP1DS,ZP1DS 100 WRITE (6,560) LDP2DS,ZP2DS 101 WRITE (6,560) LSTPMX,ZMAX 102 WRITE (6,560) LUHAG,ZHAG 103 WRITE (6,560) LPCTCC,PCTCC 104 WRITE (6,570) LGAHAS,ZTEP APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE ITERO 92 105 WRITE (6,570) LP1TRY,Z1TRY 106 WRITE (6,570) LP2TRY,Z2TRY 107 WRITE (6,570) LYPRED,N00,N00,(ZYES(I),1=1,4) 108 280 CONTINUE 109 IF (ITERF .EQ. ,2) GO TO 300 110 IF (NEQC .EQ. NDEPV) GO TO 290 11.1. WRITE (6,580) LNAC, NAC 112 WRITE (6,580) LIAC, (IAC (I) ,1=1, NAC) 11.3 290 CONTINUE 114 WRITE (6,560) LDP1DS,DP1DS 115 WRITE (6,560) LDP2DS,DP2DS 116 WRITE (6,560) LSTPMX,STPHAX 117 WRITE (6,560) LUMAG,UHAG 118 WRITE (6,560) LDUMAG,DUMAG 119 WRITE (6,570) LGAMAS,STEP 120 WRITE (6,570) LP1TRY,P1TRY 121 WRITE (6,570) LP2TRY,P2TRY 122 FOO = 0.0 123 WRITE (6,570) LYPRED,F00,F00,(YES(I),1=1,4) 124 C 125 C ABREVIATED OUTPUT FOR TRIAL STEPS BEGINS HERE 126 C 127 300 DO 320 1=1,NINDV 128 TEMP(I)=U(I)+GAMAS*DU(I) 129 320 CONTINUE 130 IF (NOHF.NE.O) GO TO 340 131 WRITE (6,560) LGAHA,GAHAS 132 WRITE (6,560) LDU,(DU(I),1=1,NINDV) 133 340 CONTINUE 134 WRITE (6,560) LWU, (TEMP (I) ,1=1,NINDV) 135 DO 360 1=1,NINDV 136 TEMP (I) =TEMP (I) /HU (I) 137 360 CONTINUE 138 WRITE (6,560) LU, (TEMP (I) ,1=1, NINDV) 139 380 CONTINUE 140 IF (NDEPV.EQ.O) GO TO 420 141 WRITE (6,560) LE, (E (I) ,1=1, NDEPV) 142 420 CONTINUE 143 IF (OPT.EQ.0.0) GO TO 440 144 WRITE (6,560) LWP1,P1 145 TEMP (1)=Pl/WOPT 146 WRITE (6,600) 147 WRITE (6,560) LP1,TEHP(1) 148 GO TO 460 149 440 CONTINUE 150 WRITE (6,600) 151 460 CONTINUE 152 WRITE (6,560) LP2,P2 153 IF (IDEB .EQ-.O) RETURN 154 IF (NOMF .EQ. 0) RETURN 155 CALL PAGER (1000) 156 RETURN 157 C APPPENDIX III - CONSTRAINED TARGETING ROUTINES 93 SUBROUTINE ITERO 158 C FORMATS 159 C 160 C 161 480 FORMAT (23H *** NOMINAL TRAJECTORY///) 162 500 FORMAT (22H *** ITERATION NUMBER , 12///) 163 560 FORMAT (A4,5 (E15.8,7X) ,E15.8/6 (7X, E15.8)/6(7X,E15.8)/6(7X#E15. 164 1X,E15.8) 165 570 FORMAT (A4, 5 (E15. 8#7X) ,E15.8) 166 580 FORMAT (A4,13X,5(I2,20X),12/6(20X,12)/6 (20X, 12)/6 (20X,I2) /20X,I2) 167 600 FORMAT (1H ) 168 620 FORMAT (///15H-*** TRIAL STEP) 169 END APPPENDIX III - CONSTRAINED .TARGETING ROUTINES SUBROUTINE MATPY 94 1 SUBROUTINE MATPY (A,B,C,L,M, N) 2 DIMENSION A(1), B(1), C(1) 4 DO 100 1=1,L 5 DO 100 J=1,N 6 JH1 = J-1 7 JM1M = JM1*H 8 II = JM1*L+I 9 C (II) = 0. , 10 DO 100 K=1,M 11 JJ .= (K-1) *L+I 12 KK = JH1M+K 13 100 C(II) = C(II)+A(JJ) *B(KK) 14 C 15 RETURN 16 END APPPENDIX III - CONSTRAINED TARGETING EOOTINES 95 SUBROUTINE MINMYS 1 SUBROUTINE HINMYS 2 C 3 C. MINMYS MINIMIZES THE OPTIMIZATION INDEX WHILE SATISFYING THE 4 C CONSTAINTS 5 C 6 C 7 C COMMON /SEARCH/ GOES HERE 8 C COMMON /SERVC/ GOES HERE 9 C 10 INTEGER SRCHM 11 C 12 C ITERATION INITIALIZATION 13 C 14 IF (SRCHM.NE.5) GO TO 20 15 SRCHM=1 16 IDAV=63 17 20 CONTINUE 18 DO 40 I = 1, NINDV 19 IF (MODEW .NE. 0) GO TO 40 20 D(I)=U(I)*SU(I) 21 PERT (I)=PERT (I) *8U (I) 22 40 CONTINUE 23 60 CONTINUE 24 CTHAT=COS(CONEPS (1) *RPD) 25 NEQC = 0 26 DO 70 I = 1, NDEPV 27 IAC(I) = I 28 IF (IDEPVR(I) .EQ-,0) NEQC = NEQC + 1 29 70 CONTINUE 30 NAC = NDEPV 31 IF (OPT .EQ. 0.0) GO TO 75 32 IF (SRCHM .NE. .4) GO TO 75 33 NETF =63 34 75 CONTINUE 35 C 36 C ITERATION LOOP 37 C 38 DO 100 NSTEP=1,100 39 CALL NOMINL 40 IF (NFLAG .NE.,0) GO TO 120 41 C 42 C GRADIENTS COMPUTED 43 C 44 CALL GRAD 45 IF (NFLAG.NE.O) GO TO 120 46 C 47 C DIRECTION OF SEARCH CALCULATED AND CONVERGENCE FLAG SET 48 C 49 ' CALL DELTU 50 IF (NFLAG.NE.O) GO TO 120 51 C APPPENDIX III - CONSTRAINED TARGETING ROUTINES 96 SUBROUTINE HINMYS 52 C FIND THE BEST STEP SIZE IH THE DIRECTION OF SEARCH 53 C 54 C 55 ITERF = 0 56 IF (NETF .EQ. 0) GO TO 80 57 IF (INTRBL.NE. 0) GO TO 80 58 ITERF = 2 59 CALL TRYIT1 60 IF (NFLAG.NE.O) GO TO 120 61 I = NTC + NEQC 62 IF (I .EQ. 0) GO TO 100 63 ITERF = 1 64 80 CONTINUE 65 CALL TRYIT2 66 IF (NFLAG.NE.O) GO TO 120 67 C 68 C UPDATE NOMINAL CONTOLS 69 C 70 CALL UPNOM 71 C 72 100 CONTINUE 73 120 CONTINUE 74 GAHAS = 0.0 75 CALL NOHINL 76 NFLAG = 0 77 RETURN 78 C . 79 END APPPENDIX III — CONSTRAINED .TARGETING ROUTINES SUBROUTINE NOMINL 97 1 SUBROUTINE NOMINL 2 C 3 c*** RUNS A NOMINAL TRAJECTORY, SAVING THE STATE VARIABLES AT THE 4 C BEGINNING OF EACH PHASE, CALLS TEST AND DISPLY. 5 C s . 6 C COMMON /SEARCH/ GOES HERE 7 C COMMON /SERVC/ GOES HERE 8 DMAG=SQRT(SP(U,U,NINDV)) 9 NOMF=63 10 C 11 CALL TRAJ 12 C 13 P1N0M=P1 14 P2NOM = O.O 15 IF (NDEPV.EQ.O) GO TO 40 16 NTC = 0 17 DO 20 1=1,NDEPV 18 ENOM (I)=E(I) 19 ITC (I) = 0.0 20 IF (IDEPVR (I) .EQ. 0) GO TO 10 21 IF (ENOM(I) .GT. 10.0) GO TO 20 22 ITC (I) =63 23 NTC = NTC + 1 24 10 CONTINUE 25 P2NOM = P2NOM • ENOM(I)**2 26 20 CONTINUE 27 40 CONTINUE 28 CALL DISPLY 29 NOMF=0 30 CALL TEST 31 RETURN 32 60 CALL ERROR (10HUNUSABLE ,10HNOHINAL ) 33 NFLAG--1 34 RETURN 35 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE PAD 98 1 SUBROUTINE PAD (A,B,HOD) 2 INTEGER TRACE 3 DIMENSION LOCREF(1) 4 DATA MIN, MAX /17,-1/ 5 DATA NAD /0/ . 6 IF (HOD .NE..0) GO TO 10 7 CALL DIGDIF(A,B,NDIF1) 8 X = A 9 IF (ABS(B) . GT. ABS (A) ) X = B 10 CALL DIGDIF (A+X,B+X,NDIF2) 11 NDIF = MINO(NDIF1,NDIF2) 12 IF (NDIF •GT..MAX) MAX = NDIF 13 IF (NDIF .EQ.,0) RETURN 14 IF (NDIF .LT. MIN) MIN = NDIF 15 RETURN 16 10 CONTINUE 17 IF (MIN .GE..4 .OH. MAX .GE. 6) GO TO 20 18 AS = A 19 A = A*1O.0**MINO(4-MIN,6-MAX) 20 GO TO 30 21 20 IF (MAX .LE. 6 .OR. MIN .LE.,4) GO TO 40 22 AS = A 23 A = A/10*0**HINO(MIN-4,MAX-6) 24 IF (B -EQ. ,0.0) GO TO 30 25 IF (ABS (A/B) .LT. , 1.0E-10) A = SIGN(B*1.0E-10, A) 26 30 CONTINUE 27 NAD = NAD + 1 28 IF (A .LT. ,1000.0 - .AND.., NAD .LT. .10) GO TO 50 29 A = AS 30 40 CONTINUE 31 NAD = 0 32 50 CONTINUE 33 WRITE (6,1) MIN, MAX 34 MIN = 17 35 MAX = -1 36 RETURN 37 1 FORMAT (7H0MIN =,13/ 38 1 7H MAX =,13/) 39 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE PAGER 99 1 SUBROUTINE PAGER (N) 2 DIMENSION HEADER(25) 3 C 4 DATA NPAGE /0 / 5 DATA KOUNT /0/ 6 DATA HEADER /25*4H / 7 C 8 K0UNT=K0UNT+N 9 IF (KOUNT.LE.60 ) RETURN 10 IF (KOUNT.GT.500) GO TO 20 11 K0UNT=N+3 12 10 NPAGE=NPAGE+1 13 WRITE (6,1000) NPAGE, (HEADER (L)#L=1,25) 14 1000 FORMAT (1H1,120X,5HPAGE ,I4,/1X,25A4/) 15 RETURN 16 20 CONTINUE 17 K0UNT=3 18 GO TO 10 19 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES 100 SUBROUTINE PGM 1 SOBROOTINE PGM 2 C COMMON /SEARCH/ GOES HERE 3 C COMMON /SERVC/ GOES HERE 4 C 5 CALL UPDATS 6 C THIS CARD SEEMS UNNECESSARY TO HE, 7 C BOT I LEAVE IT HERE JOST INCASE THERE HASD 8 C SOME REASON FOR IT.. 9 C IF (OPT .EQ. 0.0) GO TO 40 10 IF (INTRY1 .NE. 0) GO TO 40 11 IF (INTRBL .NE. 0) GO TO 40 12 CALL DROP 13 40 CONTINUE 14 IF (NAC .LE. , NINDV) GO TO 50 15 C 16 C LEAST SQUARES STEP 17 C 18 CALL BTW (SMAT,SMAT,SSTI,NINDV,NDEPV,NINDV) 19 CALL INVM(SSTI,NINDV) 20 CALL ABT (SSTI,SHAT,EHAP,NINDV,NINDV,NDEPV) 21 INTRBL = 63 22 C 23 50 CONTINUE 24 K = 0 25 DO 60 I = 1, NINDV 26 G2 (I) =0.0 27 IF (NAC .EQ. 0) GO TO 60 28 DO 55 J = 1, NAC 29 K = K • 1 30 EMAP(K) = -EMAP(K) 31 L = IAC(J) 32 EA(J) = ENOM(L) 33 M = J + (I-1)*NAC 34 G2(I) = G2(I) + FP2*EA(J) *SMAT (M) 35 55 CONTINUE 36 60 CONTINUE 37 P2NOM = 0.0 38 IF. (NAC •NE. 0) P2NOH = SP (EA,EA,NAC) 39 IF (NETF .EQ. 0) GO TO 70 40 IF (INTRBL .NE. 0) GO TO 70 41 IF (INTRY1 .NE. 0) GO TO 70 42 C 4 3 CTHA =0.0 44 IF (ABS(PGIMAG) -GT. 1.0E-8) CTHA=SP(G1, PG1,NINDV) / (G 1MAG*PG1 MAG) 45 IF (CTHA .GT. 1.0) CTHA = 1.0 46 IF (ISTART .NE. 0) GO TO 63 47 IF (NACS .NE. NAC) GO TO 62 48 IF (IDEPVR(20) . NE. ,0) GO TO 62 49 IF (NAC .EQ.,0) GO TO 65 50 DO 611 = 1, NAC 51 IF (IACS(I) .NE. IAC(I) ) GO TO 62 52 61 CONTINUE APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE PGH 101 53 GO TO 65 54 62 CONTINUE 55 ISTART = 63 56 63 CONTINUE 57 C RESTART DAVIDON 58 L = 0 59 DO .64 1 = 1 , NINDV 60 DU(I) = -PG1(I) 61 S(I) = DU(I) 62 GP (I) = PG1 (I) 63 L = L + 1 64 H(L) = FP1 65 64 CONTINUE 66 GO TO 66 67 65 CONTINUE 68 CALL DGH (PG1) 69 66 CONTINUE 70 NACS = NAC 71 DO .67 I = 1, NAC 72 IACS (!) = IAC(I) 73 67 CONTINUE 74 GO TO 100 75 C 76 70 CONTINUE 77 CALL -HATPY (EHAP,EA,DU,NINDV,NAC,1) 78 100 CONTINUE 79 RETURN 80 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES 102 SUBROUTINE PPT 1 SUBROUTINE PPT (P,C,H,N,S) 2 C 3 C*****P(MXN) *TRS(P) = C(HXM) a c 5 DIMENSION P(1), C(1) 7 DO 200 1=1,H 8 DO 200 J=I,M 9 IJ = I* (J-1) *M 10 C(IJ) = C(IJ) *S 11 DO 100 K=1,N 12 KM1SM = (K-1)*M 13 IK = I+KM1SM 14 JK = J+KM1SM 15 100 C(IJ) = C(IJ) +P(IK) *P(JK) 16 JI = J+(I-1) *M 17 200 C(JI) = C(IJ) 18 C 19 RETURN 20 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES 103 SUBROUTINE REVISE 1 SUBROUTINE REVISE 2 C COMMON /SEARCH/ GOES HERE 3 C COMMON /SERVC/ GOES HERE 4 C 5 10 CONTINUE 6 ICD = 0 7 IF (NAC .EQ.,0) GO TO 200 8 TEMP(1)= PG1MAG/G1MAG 9 IF (TEMP(1) .GT. PGEPS) GO TO 200 10 CALL MATPY (SSTI,SMAT,PROJ,NAC,NAC,NINDV) 11 CALL MATPY (PR0J,G1,TEMP(1),NAC,NINDV,1) 12 DO 60 I = 1, NAC 13 STEMP (1) = 0.0 14 DO 4 0 J = 1, NINDV 15 K = I • (J-1) *NAC 16 STEMP(1)= STEMP (1) + SMAT(K)**2 17 40 CONTINUE 18 TEMP (I) = TEMP (I) *SQRT(STEMP (1) ) 19 WRITE (6,1) I, TEMP(I) 20 60 CONTINUE 21 STEMP(1)= 0.0 22 DO 80 1 = 1 , NAC 23 J = I AC (I) 24 IF (IDEPVR(J) .EQ. 0) GO TO 80 25 IF (TEMP (I) .GT., STEMP (1)) GO TO 80 26 K = IAC(I) 27 STEMP (1)= TEMP (I) 28 80 CONTINUE 29 IF (STEMP(I).GE. -1.0E-6) GO TO 200 30 ICD = K 31 ITC(ICD) = 0 32 NTC = NTC - 1 33 CALL UPDATS 34 200 CONTINUE 35 RETURN 36 1 . FORMAT (13H-***COMPONENT,I2,3H = ,E21.14) 37 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES 104 SUBROUTINE RSEARC 1 SUBROUTINE RSEARC 2 C 3 C RSEARCH - READS NAMELIST SEARCH 4 C 5 C 6 C COHHON /SEARCH/ GOES HERE 7 C COMMON /SERVC/ GOES HERE 8 C 9 INTEGER SRCHM 10 REAL INDPH 11 C 12 NAMELIST/SEARCH/ 13 1 CONEPS#CONSEX,FITERR,GAHAX ,IDEB ,IDEPVR, KASE ,MAXITR,MODES 14 2,NDEPV #NINDV ,OPT #PCTCC ,PERT ,PGEPS , P2MIN ,SRCHM ,STMINP 15 3#U ,HCON rHOPT ,HU 16 C 17 READ (5#SEARCH) 18 RETURN 19 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE SDH SUBROUTINE SDH C C*** CALCULATES DU BY THE STEEPEST DESCENT METHOD (-G2) C C COMMON /SEARCH/ GOES HERE C COHMOH /SERVC/ GOES HERE C DO 20 1=1,NINDV DU (I)=-G2(I) 20 CONTINUE RETURN END APPPENDIX III - CONSTRAINED TARGETING .ROUTINES FUNCTION SP 106 1 FUNCTION SP (X,Y,N) 2 c * * * * * * * * * * * * * * * * * * * * * * 3 C * *** FUNCTION SP *** 4 C * PURPOSE COMPUTES THE INNER PRODUCT OF TWO N-DIMENSIONAL VECTOR 5 C * 6 C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 7 C 8 C 9 C **************************** 10 C * DIMENSION STATEMENT(S) * 11 C **************************** 12 C 13 C 14 DIMENSION 15 1 X(1) ,1(1) 16 C 17 SP=0.0 18 DO 10 1=1,N 19 SP=SP+X(I) *Y(I) 20 10 CONTINUE 21 RETURN 22 END i APPPENDIX III ~ CONSTRAINED TARGETING ROUTINES SUBROUTINE TEST 107 1 SUBROUTINE TEST 2 C 3 C THIS ROUTINE TESTS TO SEE IF THE CONVERGENCE CRITERIA HAVE BE 4 C MET 5 C 6 C COMMON /SEARCH/ GOES HERE 7 C COMMON /SERVC/ GOES HERE 8 C 9 INTEGER SRCHM 10 IF (ABS (P2NOM) .GT. P2MIN) GO TO 50 11 IF (NETF .EQ. 0) GO TO 20 12 IF (CTHA.GT.CTHAT) GO TO 60 13 C 14 C CONVERGENCE CRITERIA MET 15 C 16 20 CONTINUE 17 WRITE (6,160) 18 40 NF1AG=63 19 RETURN 20 50 CONTINUE 21 IF (IOPT .NE. 0) GO TO 60 22 IF (NETF .EQ. 0) GO TO 60 23 INTRBL = 63 24 60 IF (NSTEP.LT.3) GO TO 120 25 C 26 C COMPARE U,P1,P2,G2, STHT THE PREVIOUS ITERATION 27 C 28 RATIO(1)=ABS((UMAG-OLDU)/UMAG) 29 IF (P1NOM .NE. 0.0) RATIO (2) =ABS ( (P1NOH-OLDP1) /P1NOM) 30 RATIO(3)=ABS((P2NOM-OLDP2)/P2NOM) 31 IF (G2MAG .NE. 0) RATIO(4) = ABS((G2MAG— OLDG2)/G2MAG) 32 NVOTE=0 33 DO 80 1=1,4 34 IF (RATIO(I) .GT.CONEPS (1 + 1) ) NVOTE=1 35 80 CONTINUE 36 IF (KREEP.NE.O.OR.NVOTE.NE.0) GO TO 100 37 C 38 C IF THERE IS NO CHANGE 2 ITERATIONS IN A ROW, STOP 39 C 40 WRITE (6,180) 41 GO TO 40 42 100 CONTINUE 43 KREEP=NVOTE 44 120 OLDU=UMAG 45 0LDP1=P1N0M 46 OLDP2=P2NOM 47 OLDG2=G2MAG 48 IF (NSTEP.LE.MAXITR) GO TO 140 APPPENDIX III - CQNSTRAINED TARGETING ROUTINES SUBROUTINE TEST 108 49 50 51 140 52 53 C 54 160 55 180 56 57 200 58 NFLAG=-1 WRITE (6,200) CONTINUE RETURN FORMAT (20H-*** FORMAT (55H-*** D FORMAT (20H-*** END PROBLEM SOLVED ) NO CHANGE IN STATE DURING 2 CONSECUTIVE ITERATIONS ITERATION LIMIT) APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE THPH 109 1 SUBROUTINE THPH(X,Y,XMIN,YMIN) 2 DIMENSION X(3), Y (3) 3 X12 = X(1) - X(2) H X23 = X(2) - X(3) 5 X31 = X(3) - X(1) 6 ANUH = X23*Y(1) + X31*Y(2) + X12*Y(3) 7 XS12 = X12*(X(1)+X(2)) 8 XS23 = X23* (X (2) *X (3) ) 9 XS31 = X31*(X(3)+X(1)) 10 DENOM = XS12*X(3) + XS23*X(1) + XS31*X(2) 11 A = -ANUM/DENOH 12 IF (A .GT. 0.0) GO TO 10 13 XHIN = 1..0B10 14 YHIN = -1.0E10 15 RETURN 16 10 BNUM = XS23*Y(1) • XS31*Y(2) + XS12*Y(3) 17 B = BNUH/DENOH 18 XMIN = -.5*B/A 19 C = Y(1) - X(1)*(B + X(1)*A) 20 YHIN = C • XMIN*(B • XHIN*A) 21 RETURN 22 END APPPENDIX III - CONSTRAINED TARGETING ROOTINES SUBROUTINE THPOSM 110 1 SUBROUTINE THPOSM(X,Y,DYDX1,XMIN,¥HIN,IEBR) 2 C 3 C INPUTS 4 C X=ARRAY OP ABSCISSA VALUES OF THREE SAMPLE POINTS 5 C Y=ARRAY OF ORDINATE VALUES OF 'THREE SAMPLE POINTS 6 C DYDX1=SLOPE AT FIRST SAMPLE POINT 7 C OUTPUTS 8 C XMIN=ABSCISSA VALUE OF MINIMUM OF CUBIC FITTED THROUGH 9 C THREE SAMPLE POINTS AND THE SLOPE AT FIRST POINT 10 C YMIN=ORDINATE VALUE OF MINIMUM OF ABOVE CUBIC 11 C IERR=FLAG WHOSE NON-ZERO VALUE INDICATES THAT TWO OF THE 12 C GIVEN X VALUES ARE IDENTICAL. 13 C 14 DIMENSION 15 1 X(3) ,Y(3) ,A(4) 16 C 17 DATA 18 1 ONE /1. / 19 2 ,TWO /2. / 20 3 ,THREE /3. / 21 C 22 DATA 23 1 NO /0 / 24 2 ,N1 /1 / 25 C 26 IERR=N0 27 IF (X (3) —X (2) ) 30,80,40 28 30 CONTINUE 29 AMBDA=X(2)-X(1) 30 IF (AMBDA) 70,80,70 31 80 CONTINUE 32 IERR=N1 33 GO TO 60 34 70 CONTINUE 35 ALPHA=(X(3)-X(1) )/AMBDA 36 FALAM=Y(3) 37 FLAM=Y(2) 38 GO TO 50 39 40 CONTINUE 40 AMBDA=X<3)-X(1) 41 IF(AMBDA) 90,80,90 42 90 CONTINUE 43 ALPHA=(X(2)-X(1))/AMBDA 44 FALAM=Y (2) 45 FLAM=Y (3) 46 50 CONTINUE 47 IF (ALPHA) 100,80,100 APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE THPOSM 111 100 CONTINUE F0=Y(1) A(1)=F0 F0P=DYDX1 A{2)=F0P ALPHAS=ALPHA**2 APSFLH=ALPHAS*FLAM OPALP=ONE+ALPHA OHALP=ONE-ALPHA APSAHS=AMBDA**2*ALPHAS ALFOP=ALPHA*AHBDA*F0P A{«)=((APSFLM-FALAM)/OMALP+ALF0P+ F0*OPALP)/(APSAMS*AMBDA) & (3) = ( (FALAH-ALPHA*APSFLM) /OHALP-ALF0P*OPALP-F0*(OPALP+ALPHAS 1 /APSAHS CALL CDBHIN (A,XHIN,YHIN) 60 CONTINUE RETURN END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SOBROOTINE TPOSM 112 1. SUBROUTINE TPOSM(X,Y,D¥DX1,XMIN,YMIN) 2 DIMENSION X (2) , Y (2) 3 X2MX1 = X(2) — X(1) 4 A = ( ¥ ( 2 H Y(1) )/X2MX1**2 - DYDX1/X2MX1 5 IF (A .GT. ,0.0) GO TO 10 6 XMIN = 1.0E10 7 YMIN = -1.0E10 8 RETURN 9 10 A2 = 2.0*A 10 B = DYDX1 - A2*X (1) 11 XMIN = -B/A2 12 C = Y(1) - X(1)*(B • X(1)*A) 13 YMIN = C • XMIN*(B + XMIN*A) 14 RETURN 15 END APPPENDIX I I I — CONSTRAINED TARGETING ROUTINES 113 SUBROUTINE TRAJ 1 p SUBROUTINE TRAJ 3 C*** RUNS THE TRAJECOTRY PHASE BY PHASE, SETING • IN THE INDEPENDENT 4 c AND SAVING THE DEPENDENT VARIABLES AT THE APPROPRIATE PHASES. 5 c THE CP TIHE FOR THE TRAJECTORY AND OTHER OPTIMIZATION DATA ARE OUT 6 c PUT 7 c 8 c COHHON /SEARCH/ GOES HERE 9 c COMMON /SERVC/ GOES HERE 10 c 11 INTEGER SRCHM 12 COMMON /COM/ C(12) 13 DIMENSION P(13), PO (13) , PD(13), Q(13) 14 DIMENSION DEPTL (25) 15 c 16 DATA DEPTL /25*.01/ , 17 DATA PO /O.ODO, 18 1 5.0D1,7.1D1,1.12D2,3.0D0, 19 2 5.0Dl,7.iD1,1.12D2,3.0D0, 20 2 2.6D1,4.5D1,9.9D1,1.5D1/ 21 CALL TIME (1,0,1) 22 TEMPT = I 23 DO 999 1 = 1 , NINDV 24 TEMP (I) = (U(I) + GAMAS*DU(I))/HO(I) 25 999 CONTINUE 26 C GO TO (1010,1020,1030), KASE 27 C1010 CONTINUE 28 C E(1) = (TEMP(3)-FPP5*(TEHP(1)**2+ 29 FP2*TEMP(2) **2) )/DEPTL (1) C E(2) = (TEMP (2) - FP1)/DEPTL(2) 30 C E(3) = (TEMP (1) - FPP5)/DEPTL (3) 31 C E(4) = (TEMP(1) + TEMP(2) - FP2) /DEPTL (4) 32 C E(5) = (FP4- (TEMP(1)+TEMP (2)) )/DEPTL (5) 33 C P1 = TEMP (3) 34 C GO TO 2000 35 C1020 CONTINUE 36 C E(1) = (TEMP (3) -FPP5* (TEHP(1),**2+ 37 TEMP(2) **2) )/DEPTL (1) C E(2) = (TEMP(2)-FP3+FPP5*(FP5*TEMP(1) TEMP(1) **2) )/DEPTL (2) 38 C P1 = TEMP (3) 39 c GO TO 2000 40 C1030 CONTINUE 41 C TEMP(10) = (TEHP(1)-FP1) **2 42 C TEMP(11) = (TEHP(1) +FP1)**2 43 C E(1) = (TEMP (10)-TEMP (2))/DEPTL (1) 44 C . E(2) = (TEMP (10) +TEMP(2))/DEPTL(2) 45 c E(3) = (TEMP (11)-TEMP (2))/DEPTL (3) 46 c E(4) = (TEMP(11)+TEMP(2))/DEPTL(4) 47 c E(5) = (FP1-(TEMP(1) **2 + APPPENDIX III - CONSTRAINED TARGETING ROUTINES 114 SOBROOTINE TRAJ TEHP (2) **2) ) /DEPTL (5) 48 C P1 = TEHP (1) +TEMP (2) 49 C2000 CONTINUE 50 DO 10 I = 1, 13 51 P(I) = P0(I) 52 10 CONTINOE 53 DO 20 I = 1, NINDV 54 C(I) = TEHP (I) 55 20 CONTINUE 56 PD(1) •= 1.0 57 CALL AUXRK(P,PD) 58 IF (NOMF .NE. 0) WRITE (7,1) (P (J) , J=2,13) , (PD (J) ,J=2,13) 59 DO 40 I = 1, 80 60 CALL RK(P,PD,Q,1.0E-2,13,100) 61 IF (I .NE. 1) GO TO 30 62 E(6) = -PD (2)* 1000.0 63 E(7) = -PD (3)* 1000.0 64 E(8) = -PD(4) *1000.0 65 E(9) = -PD(5)*1000.0 66 30 CONTINUE 67 IF (NOHF .NE.,0) WRITE (7,1) (P (J) , J=2,13) , (PD(J) ,J=2,13) 68 CC IF (P(J) .LT. ,1) P(J) = 0.0D0 69 C30 CONTINUE 70 CC WRITE (6,123) P 71 C123 FORMAT (1X,13D10.3) 72 40 CONTINUE 73 E(1) = P(2) - 3.OE0 74 E(2) = P(3) - 3.6E1 75 E(3) = P(4) - 1.1E1 76 E(4) = P(5) - 7.0E0 77 E(5) = AMAX1 (P (6) ,0.0)+AMAX1 (P (7) ,0.0) + AMAX1 (P(8) ,0.0) - 1.5E0 78 C E(6) = P(7) - O.OEO 79 C E(7) = P(8) - O.OEO 80 C E(8) = P(9) - O.OEO 81 E(10) = -PD(5) *1000.0 82 E(11) = -PD (3) *1000. 0 83 E(12) = -PD(4) *1000.0 84 E(13) = AMAX1 (P(10) ,0.0)+AMAX1 (P (11) ,0.0) + AMAX1 (P (12) ,0.0) - 1.5E0 85 E(14) = -PD(5) *1000.0 86 P2 = 0.0 87 IF (NAC .EQ.,0) GO TO 88 88 DO 85 I = 1, NAC 89 J = IAC(I) 90 EA(I) = E(J) 91 85 CONTINOE 92 P2 = SP (EA,EA,NAC) 93 88 CONTINUE 94 IF (SRCHM.EQ.1) P2=P1+WC0N*P2 95 IF >(NOMF .NE. ,0) GO TO 90 96 IF (IDEB.EQ.O) GO TO 100 APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE TRAJ 115 97 CALL TIME (1,0,1) 98 TEMP(1) = I -TEMPT 99 WRITE (6,120) TEMP(1) 100 CALL PAGER (N02) 101 IF (GAMAS .EQ. 1.0) GO TO 100 102 90 CONTINUE 103 CALL ITERO 104 100 CONTINUE 105 RETURN 106 C 107 1 FORMAT (12(1X,El0. 3) ) 108 120 FORMAT (25H *** TRAJECTORY CP TIME =,F7.3/) 109 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE TRYIT1 116 1 SUBROUTINE TRYIT1 2 C 3 C*** MINIMIZES NET PERFORMATNCE 4 C 5 C COMMON /SEARCH/ GOES HERE 6 C COMMON /SERVC/ GOES HERE 7 DIMENSION EPS(4) 8 C 9 INTRY1 = 63 10 IF (IOPT .EQ. 0) GO TO 40 11 C 12 C ADJUST PCTCC BASED UPON PREVIOUS ITERATION 13 C 14 IF (IHADIT.NE.0) GO TO 20 15 PCTCC = FP2*PCTOLD 16 GO TO 40 17 20 CONTINUE 18 PCTCC=PCTOLD*FPP5 19 C 20 40 CONTINUE 21 IHADIT =0 22 DP1DS=SP(DO,G1,NINDV) 23 DP2DS=SP(DU,G2,NINDV) 24 EPS(1) = STMINP 25 EPS(2) = AMIN1 (STPMAX,PCTCC*UMAG) 26 EPS (3) = CONSEX 27 EPS (4) = FITERR 28 C 29 C INITIALIZE TRIAL STEP ARRAY 30 C 31 DO 80 1=1,6 32 STEP(I)=0.0 33 PITRY(I) =0.0 34 P2TRY(I)=0.0 35 IF (I.GT. 4) GO TO 80 36 YES (I) = 0.0 37 80 CONTINUE 38 P1TRY(1)=P1N0M 39 P2TRY (1) =P2NOM 40 IF (NAC .EQ. 0) GO TO 95 41 CALL MATPY (EMAP,EA,TEMP (1) ,NINDV,NAC,1) 42 P1TRY(1) =P1TRY (1) +SP (G1,TEMP (1) , NINDV) 43 95 CONTINOE 44 STEP (2) = FPP5*EPS(2) 45 C 46 C MINIMIZE PERFORMANCE 47 C 48 CALL GENMIN (STEP,P1TRY,DP1DS,EPS,YES,MIN) 49 IF (MIN.NE.O) GO TO 100 50 WRITE (6, 180) 51 IF (ISTART .NE. 0) NFLAG = -1 52 1 = 1 53 GAMAS = 0 APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE TRYIT1 117 54 GO TO 110 55 100 CONTINOE 56 BIN = IABS(MIN) 57 GAMAS = STEP (BIN) 58 . IF (GAMAS .NE* ,EPS (2)) ISTART = 0 59 GAMASS = GAMAS 60 C 61 C SAVE RESULTS TO BE OUTPUT IN THE ITERATION SUMMARY 62 C 63 110 CONTINUE 64 PCTOLD = PCTCC*FPP5 65 IF (GAMAS .NE. STPMAX) PCTOLD = GAMAS/UMAG 66 ZP1DS=DP1DS 67 ZP2DS=DP2DS 68 ZMAG = UMAG 69 ZMAX = STPMAX 70 DO 120 J = 1,6 71 ZTEP(J) •= STEP(J) 72 Z1TRY(J) = P1TRY(J) 73 Z2TRY(J) = P2TRY(J) 74 IF (J .GT. 4) GO TO 120 75 ZYES(J) = YES(J) 76 120 CONTINOE 77 JAC(1) = NAC 78 DO 130 I = 1, NAC 79 JAC(I+1) = IAC (I) 80 130 CONTINUE 81 C 82 C UPDATE THE CONTROLS 83 C' 84 I = MIN - 1 85 P1NOM=SAVIT(1,I) 86 P2NOM=SAVIT (2,1) 87 NTC = 0 88 DO 140 J=1,NDEPV 89 ENOM (J)=SAVIT (J+2,I) 90 IF (IDEPVR(J) .EQ. ,0) GO TO 140 91 IF (ENOH(J) .GT. ,1.0) GO TO 140 92 ITC(J) = 63 93 NTC = NTC • 1 94 140 CONTINUE 95 CALL UPNOH 96 C 97 C CALCULATE A TARGETING DU 98 C 99 IF (NFLAG .EQ. 0) CALL DELTU 100 160 CONTINUE 101 IOPT = 63 102 INTRY1 =0 103 RETURN 104 C 105 180 FORMAT (20H OPTIMIZATION UPHILL) 106 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE TRYIT2 11.8 1 SUBROUTINE TRYIT2 2 C 3 C*** MINIMIZES THE ERROR 4 C 5 C COMMON /SEARCH/ GOES HERE 6 C COMMON /SERVC/ GOES HERE 7 C 8 INTEGER SRCHM 9 DIMENSION EPS (4) 10 DP1DS=SP(DU,G1,NINDV) 11 DP2DS=SP(DU,G2,NINDV) 12 IF (INTRBL .NE. 0) GO TO 10 13 IF (NETF .EQ. 0) GO TO 10 14 DP2DS = 0.0 15 10 CONTINUE 16 INTRBL = 0 17 C 18 C INITIALIZE TRIAL STEP ARRAY 19 C 20 DO 20 1=1,6 21 P1TRY(I) =0.0 22 P2TRY(I) =0.0 23 STEP (I) =0.0 24 20 CONTINUE 25 P1TRY(1) =P1N0M 26 P2TRY (1) =P2NOM 27 EPS(1) = STMINP 28 EPS (2) = AMIN1 (STPHAX,GAHAX*DUMAG) 29 EPS(3) = CONSEX 30 EPS (4) = FITERR 31 STEP(2)=ABS(FP2*P2NOM/DP2DS) 32 C IF {SRCHM. EQ. 4) STEP (2) =DUMAG 33 IF (SRCHM. EQ-4) READ (5,1) STEP (2) 34 1 FORMAT (E10.0) 35 STEP(2) = AMIN1 (STEP(2) ,EPS(2)) 36 IF (DP2DS .EQ. .0.0) STEP (3) = AMIN1 (1.25*DUMAG,EPS(2)) 37 IF (SRCHH.NE.4) GO TO 40 38 CALL FGAMA (N02) 39 IF (P2.GE.P2HIN) GO TO 40 40 MIN=2 41 GO TO 60 42 40 CONTINUE 43 C 44 C , MINIMIZE ERROR 45 C . 46 CALL GENHIN (STEP,P2TRY,DP2DS,EPS,YES,MIN) 47 IF (MIN.NE.O) GO TO 60 48 WRITE (6,100) 49 GAMAS = 0.0 50 IF (P2NOH .LT..P2HIN) GO TO 80 51 IF (ITERF .NE. .0) GO TO 80 52 NFLAG=-1 APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE TRYIT2 119 53 RETURN 54 60 CONTINUE 55 I=IABS(MIN) 56 GAMAS=STEP(I) 57 P2=P2TRX(I) 58 IF (P2.LT.P2HIN) RETURN 59 80 CONTINUE 60 INTRBL=63 61 IH6DIT=1 62 RETURN 63 C 64 100 FORMAT (27H DIRECTION OF SEARCH UPHILL) 65 END APPPENDIX III - CONSTRAINED TARGETING .ROUTINES SUBROUTINE UNITDU 120 1 SUBROUTINE UNITDU 2 C 3 C UNITIZE DU 4 C 5 C COMMON /SEARCH/ GOES HERE 6 C COMMON /SERVC/ GOES HERE 7 C 8 DUMAG=SQRT(SP(DU,DU,NINDV) ) 9 IF (DUHAG .EQ. ,0.0) RETURN 10 DO 20 1=1,NINDV 11 DD (I)=DU(I) /DUMAG 12 20 CONTINUE 13 RETURN 14 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE UPDATS 121 1 SUBROUTINE UPDATS 2 C COMMON /SEARCH/ GOES HERE 3 C . COMHON /SERVC/ GOES HERE 4 C 5 H = 0 6 NAC = NEQC + NTC 7 DO 60 1 = 1 , NDEPV 8 IF (IDEPVR (I) . EQ. 0) GO TO 40 9 IF (ITC(I) .EQ. 0) GO TO 60 10 40 CONTINUE 11 M = M + 1 12 IAC(H) = I 13 DO 50 J = 1 ,NINDV 14 K = M • ( J - 1 )*NAC 15 L = I + ( J - 1 )*NDEPV 16 SHAT(K) = ACOB(L) 17 50 CONTINUE 18 60- CONTINUE 19 C 20 C CALCULATE THE PROJECTED GRADIENT 21 C 22 IF (NAC .NE. ,0) GO TO 70 23 DO 65 I = 1, NINDV 24 PG1 (I) = 61.(1)-25 65 CONTINUE 26 GO TO 120 27 70 CONTINUE 28 CALL PPT (SMAT,SSTI,NAC,NINDV,0.0) 29 CALL INVH (SSTI,NAC) 30 CALL BTB (SMAT,SSTI,EHAP,NINDV,NAC,NAC) 31 IF (OPT .EQ. ,0.0) RETURN 32 CALL MATPY . (EHAP,SMAT,PROJ,NINDV,NAC,NINDV) 33 IF (NAC .LE. NINDV) GO TO 75 34 CALL ZERO (PG1,NINDV,N01) 35 GO TO 120 36 75 CONTINUE 37 DO 100 I = 1, NINDV 38 PG1 (I) = 0.0 39 DO 80 J = 1, NINDV 40 K = I + (J-1) *NINDV 41 TEMP (1)= 0.0 42 IF (I .EQ. J) TEMP(1)= 1.0 43 PG1 (I) = PG1 (I) + (TEHP (1).-PROJ (K) ) *G1 (J) 44 80 CONTINUE 45 100 CONTINUE 46 120 CONTINUE 47 PG1MAG = SQRT(SP(PG1,PG1,NINDV)) 48 RETURN 49 END APPPENDIX III - CONSTRAINED TARGETING ROUTINES SUBROUTINE UPNOM 122 1 SUBROUTINE UPNOM 2 C 3 C COMMON /SEARCH/ GOES HERE H C COMMON /SERVC/ GOES HERE 5 C 6 DO 20 I=1rNINDV 7 U(I)=U(I)+DU(I)*GAHAS 8 DU(I)=0.0 9 20 CONTINUE 10 IF (IDAV .NE. 0) GAMASS = GAMAS/DUMAG 11 . GAMAS=0 12 RETURN 13 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES 123 SUBROUTINE WUCAL 1 SUBROUTINE WUCAL 2 C 3 C*** CALCULATES THE WEIGHTING FOR THE CONTROLS 4 C 5 C COMMON /SEARCH/ GOES HERE 6 C COMMON /SERVC/ GOES HERE 7 INTEGER SRCHM 8 C 9 IF (NSTEP.NE.1) RETURN 10 IF (INTRY1 .NE. 0) RETURN 11 IF (NINDV.EQ.0) RETURN 12 IF (MODEW.EQ.O) RETURN 13 G2MAG=SQRT(SP(G2,G2,NINDV) ) 14 C 15 DO 200 1=1,NINDV 16 IF (NEQC .EQ..0) GO TO 20 17 GO TO (20,40,80,120,140,150) , MODEW 18 C . 19 C UNITIZED CONTROL WEIGHTING 20 C 21. 20 CONTINUE 22 WU (I) = 1.0 23 IF (U(I).EQ.0.0) GO TO 160 24 WU (I) =FP1/ABS (U (I) ) 25 GO TO 160 26 C 27 C DICKHAN WEIGHTING (NOT RECOMMENDED) 28 C 29 40 CONTINUE 30 TEMP (1)= 0.0 31 DO 60 J=1,NDEPV 32 IR1-J+(1-1) *NDEPV 33 TEMP(1)= TEMP(1)+ ABS (ACOB (IRI) ) 34 60 CONTINUE 35 WU (I) = TEMP (1) *WU (I) 36 GO TO 160 37 C STEINHOFF WEIGHTING (DEFINITELY NOT ' RECOMMENDED) 38 C 39 80 CONTINUE 40 WU (I) = 1.0 41 IF (U(I) .EQ.0-0) GO TO 160 42 WU (I) =0.0 43 DO 100 J=1,NDEPV 44 IR1=J+(1-1) *NDEPV 45 WU(I)=WU(I)+ABS((ACOB(IR1)*E(J))/U(I)) 46 100 CONTINUE 47 GO TO 160 48 C 49 C GRADIENT WEIGHTING 50 C 51 120 - CONTINUE 52 WU(I)=ABS(G2(I))/G2HAG APPPENDIX III - CONSTRAINED TARGETING ROUTINES 124 SUBROUTINE WUCAL r. 53 GO TO 160 54 C 55 C AVERAGED GRADIENT AND CONTROL WEIGHTING 56 C 57 140 CONTINUE 58 IF (62.(1) .EQ.0-0) GO TO 20 59 TEHP (1)=FP1/G2 (I) 60 WU(I) = (FP10*U(I)+.1*TEMP(1))/(U(I)**2+ TEHP(1)**2) 61 GO TO 160 62 C 63 C HOUSTON WEIGHTING 64 C 65 150 CONTINUE 66 TEMP (1)= 0.0 67 DO 155 J ..= 1, NDEPV 68 K = J. • (1-1) *NDEPV 69 TEHP(1) = TEMP(1) + ALOG10 (ACOB (K) **2/ACOB (J) **2 70 155 CONTINUE 71 WU(I) = SQHT (FP10**TEMP(1)/WU (I)) 72 C 73 C WEIGHT THE CONTROL,ITS PERTURBATION, ITS LIMITS, AND ITS GRADIENTS 74 C 75 160 CONTINUE 76 TEMP (1)=WU (I) 77 U(I)=U(I) *TEMP(1) 78 PERT (I) ==PERT (I) *TEMP (1) 79 G1 (I)=G1 (I)/TEHP(1) 80 G2(I)=G2(I)/TEHP(1) 81 PG1 (I) =PG1 (I)/TEMP (1) 82 DO 180 J=1,NDEPV 83 IR1=J+(1-1) *NDEPV 84 ACOB(IR1) = ACOB (IR1)/TEHP(1) 85 180 CONTINUE 86 200 CONTINUE 87 UHAG = SQRT(SP(U,U,NINDV)) 88 RETURN 89 END APPPENDIX III — CONSTRAINED TARGETING ROUTINES SUBROUTINE ZERO 125 1 SUBROUTINE ZERO (A,H,N) 2 DIMENSION A(1) 4 MN = M*N 5 DO 100 1=1,MN 6 100 A (I) = 0. , 7 C 8 RETURN 9 END APPPENDIX IV - CATERPILLAR MODEL CODE MAIN PROGRAM 126 1 C ******************************** 2 C PROGRAM CATS 3 Q ******************************** 4 C 5 C OHIT 0 HAS THE GENERAL INPUTS 6 C . UNIT 1 HAS THE WEATHER DATA IN THE FORM 7 C OF IYEAR,IDAY,IWEATHER 8 C UNIT 2 HAS THE MICRO CLIMATE IN FORM 6012 9 C UNIT 3 HAS THE TREE HEIGHTS IN FORM 6012 10 C UNIT 4 HAS THE BOUNDRIES OF THE INHABITABLE 11 C AREAS ON THE GRID 12 C UNIT 6 IS THE OUTPUT FILE FOR PRINTING 13 C UNIT 7 IS THE OUTLINE COORDINATES FOR PLOT-14 C TING 15 C UNIT 8 IS THE OUTPUT FILE FOR PLOTS 16 C 17 C 18 C 19 C FOR THE ITH, JTH GRID SQUARE: 20 C 21 C G(1,I,J) IS THE NUHBER OF TYPE I CATS 22 C G(2,I,J) IS THE NUMBER OF TYPE II CATS 23 C G(3,I,J) IS THE NUHBER OF TYPE III CATS 24 C G(4,I,J) IS THE NUMBER OF TYPE IV CATS 25 C G(5,I,J) IS THE TOTAL NUMBER OF CATS 26 C G(6,I,J) IS THE DERIVATIVE OF TYPE I 27 C. G(7,I,J) IS THE DERIVATIVE OF TYPE II 28 C G(8,I,J) IS THE DERIVATIVE OF TYPE III 29 C G(9,I,J) IS THE DERIVATIVE OF TYPE IV 30 C G(10,I,J) IS THE MICRO CLIMATE NUMBER 31 C G(11,I,J) IS THE TREE HEIGHT NUMBER 32 C G(12,I,J) IS THE I OF PREVIOUS SQUARE 33 C 6 (13,1,J) IS THE J OF PREVIOUS SQUARE 34 C G(14,I,J) IS THE I OF FOLLOWING SQUARE 35 C G(15,I,J) IS THE J OF FOLLOWING SQUARE 36 C 37 C 38 C 39 C 40 INTEGER*2 IBND 41 COMMON /CM/ 42 1 IDAY ,IYEAR , NDAY 43 2 ,NYEAR ,IHEAD ,JHEAD 44 3 ,IMAX ,JMAX ,IODAY 45 4 ,IOPLOT ,IOSAVE ,IOTOTL 46 5 ,IOTYPE ,IWPT ,MIX 47 6 , MAX ,MIY , MAY 48 7 ,ITURN ,ITILT ,MAXPOP 49 8 ,DX #DY ,SCALE 50 9 , PHAX ,C1 (4) ,C2(4) 51 0 ,C3(4) ,C4 (4) , DTAB (3) 52 A ,ETAB(5,3) ,CTAB(28) , STAB (28) 53 B ,NBD(110) ,XB (200) ,YB(200) APPPENDIX IV - CATERPILLAR MODEL CODE MAIN PROGRAM 127 54 C ,IiETH (3,500) ,IBND(2,1400) ,G (15,60, 110) 55 DIMENSION IG (15,60,110) 56 EQUIVALENCE (IG,G) 57 C 58 CALL DFAOLT (»0=TOIT:MAT.SAN.O •) 59 CALL DFAULT (*1=TOIT:MAT.REALH.1 ») 60 CALL DFAOLT (•2=TOIT:MAT.MICRO.2 •) 61 CALL DFAOLT (•3=T0IT:MAT.TREEH.3 •) 62 CALL DFAOLT (•4=T0IT:MAT.0LINE.4 •) 63 CALL DFAOLT (*6=*SINK* •) 64 CALL DFAOLT (•7=T0IT:MAT.TRACED ») 65 CALL DFAOLT (»8=-CAT.PLOTS •) 66 CALL READ 67 IF (IOPLOT .NE..0) CALL PLOTI 68 ITEAR = 1 69 IHPT = 1 70 20 CONTINOE 71 CALL INITAL 72 IF (IODAT .LE.,0) CALL OOTPUT 73 IF (ITEAR .GT. NYEAR) GO TO 100 74 40 CONTINOE 75 CALL STEP 76 IF (IDA! .EQ.IABS(IODAT) .08.. IODAY .EQ. 0) CALL OOTPOT 77 IDAT = IDAT + NDAT 78 IF (IDAT ..LE. 80) GO TO 40 79 IF (IODAY .LT. 0) CALL OOTPOT 80 CALL DSPRS 81 IYEAR = IYEAR • 1 82 GO TO 20 83 100 CONTINOE ' 84 STOP 85 END APPPENDIX IV - CATERPILLAR 130DEL CODE BLOCK DATA 128 1 BLOCK DATA 2 C 3 C GENERAL CONSTANT INITIAL VALOES 4 C 5 C DTAB IS THE DISTANCE ADULT MOTHS CAN FLY 6 C ETAB IS THE NUMBER OF EGGS AN ADULT LAYS 7 C STAB ARE SINES FOR 0 TO 90 DEGREES 8 C CTAB ARE COSINES FOR 0 TO 90 DEGREES 9 C 10 C COMMON /CM/ GOES HERE 11 DATA C1 ./.00996,. 120,.307,. 00935/ 12 DATA C2 /-.000685,-.000126,-.00055,-.000542/ 13 DATA C3 /-.000151,-.P0385,-.0100,-.000289/ 14 DATA C4 /.01>.01,-01,.01/ 15 DATA DTAB / 20.0, 12.0, 5.0 / 16 DATA ETAB / 236.0,50.6,71.6,112.0,3.0,181.0, 21.0,46.0,101.0,13.0, 17 1 140.0,7.0,18.0,85.0,30.0/ 18 DATA STAB / 1.0OO000O0E-10, .58144829E-19 01, .1 609291E+00, .17364818E+00, .23061587E+ 20 00, .28680323E+00, * .34202014E+00, .39607977E+ 21 00, .44879918E+00, * .5O000000E+00, .54950898E+ 22 00, .59715859E+00, * .64278761E+00, .68624164E+ 23 00, .72737364E+00, .76604444E+00, .80212319E+ 24 00, .83548781E+00, * .86602540E+00, .89363264E+ 25 00, .91821611E+00, * .93969262E+00, .95798951E+ 26 00, .97304487E+00, * .98480775E+00, .99323836E+ 00, .99830816E+00, 27 .10000000E+01 / 28 DATA CTAB / .10000000E+01, .99830816E+ 00, .99323836E+00, 29 * .98480775E+00, .97304487E+ 00, .95798951E+00, 30 * .93969262E+00, .91821611E+ 00, .89363264E+00, 31 * .86602540E+00, .83548781E+ 00, .80212319E+00, 32 * .76604444E+00, .72737364E+ 00, .68624164E+00, 33 * .64278761E+00, .59715859E+ 00, .54950898E+00, 34 * .50000000E+00, .44879918E+ 00, .39607977E+00, 35 * .34202014E+00, .28680323E+ APPPENDIX I? - CATERPILLAR MODEL CODE BLOCK DATA 129 00, .23061587E+00, 36 * . "J7364818E+00, . 11609291E+ 00, .58144829E-01, 37 * 1.00000000E-10 / 38 END APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE READ 130 1 SUBROUTINE READ 2 C 3 C READS THE INPUT TO THE PROGRAM 4 C 5 C COMMON /CM/ GOES HERE 6 DIMENSION IMP1 (60) , IMP2 (60) 7 C IMAX MAXIMUM GRID SIZE IN Y 8 C JHAX MAXIMUM GRID SIZE IN X 9 C IDAY CURRENT DAY OF THE YEAR 10 C (INITIALLY, THE DAY ON WHICH 11 C TO START THE SIMULATION) 12 C NDAY NUHBER OF DAYS AT A TIME TO 13 C BE INTEGRATED 14 C NYEAR NUMBER OF YEARS TO SIMULATE 15 C (SIMULATION ENDS WHEN IYEAR 16 C ] NYEAR OR THE FIRST DAY OF 17 C NYEAR + 1) 18 C IODAY ] 0 THE DAY TO OUTPUT 19 C =0 OUTPUT EVERY DAY 20 C [ 0 OUTPUT FIRST AND LAST DAY AND 21 C IABS OF IODAY 22 C IOPLOT = 2 CALCOMP PLOT LINED 23 C =1 CALCOMP PLOT OUTLINED 24 C =0 PRINTER DIAGRAM 25 C =-1 ADAGE PLOT OUTLINED 26 C =-2 ADAGE PLOT LINED 27 C IOSAVE = 0 PLOT IMMEDIATELY 28 C =0 SAVE IT PLOTS ON FILE 8 29 C IOTOTL = 0 DON'T OUTPUT TOTALS AS NUMBERS 30 C =0 OUTPUT TOTALS AS NUMBERS 31 C IOTYPE = 0 NO TYPES OUTPUT 32 C =1 TYPE 1 33 C =2 TYPE 2 34 C =3 TYPE 3 35 C =4 TYPE 4 36 C =5 TOTAL OF THE TYPES 37 C MIX MIN X AFTER ROTATION 38 C MAX MAX X AFTER ROTATION 39 C MIY MIN Y AFTER ROTATION 40 C MAY MAX Y AFTER ROTATION 41 C MAXPOP MAX POP ON ANY GRID SQUARE 42 C (USED FOR SCALING PLOTS) 43 C ITORN ANGLE RIGHT OF STRAIGHT TO 44 C VIEW PLOTS 45 C ITILT ANGLE ABOVE HORIZONTAL TO 46 C VIEW PLOTS 47 CALL FREAD (0,•171:•,IMAX,JMAX,IDAY,NDAY, NYEAR,IODAY,IOPLOT,IOSAVE 48 1 . ,IOTOTL,IOTYPE,MIX,MAX,MIY,MAY, MAXPOP,ITURN,ITILT) 49 CALL ZERO 50 CALL FREAD (-2,*ENDFILE',1) 51 5 CONTINUE APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE READ 131 52 CALL FREAD (0,•• 21: • ,1, J,S20) 53 CALL FREAD (•*»,•REAL VECTOR:1,G(1,1,3),4) 54 G(5,I,J) = 0 55 DO 10 K = 1, 4 56 G(5,I,J) = G(5,I,J) + G(K,I,J) 57 10 CONTINUE 58 GO TO 5 59 20 CONTINUE 60 1 = 1 61 30 CONTINUE 62 CALL FREAD (1,'INTEGER VECTOR:f,IHETH(1,I), 3,S60) 63 I = I + 1 64 GO TO 30 65 60 CONTINUE 66 DO .100 I = 1, IMAX 67 READ (2,1) (IMP1 (J) ,J=1,JMAX) 68 READ (3,1) (IMP2 (J) , J=1, JHAX) 69 DO 80 J = 1, JMAX 70 IF . (IHP2 (J) .EQ. .0) IMP2 (J) = -2 71 IF (IHP2(J) .EQ. 4) IHP2 (J) =-2 72 IF (IMP2 (J) . EQ. , 3) IMP2 (J) = 0 73 IMP2(J) = IMP2(J) • 1 74 IG(10,J,I) = IMP1 (J) 75 IG(11,J,I)= IMP2(J) 76 80 CONTINUE 77 100 CONTINUE 78 K = 0 79 DO 200 1 = 1 , IHAX 80 READ (4,2) N 81 . NBD(I) = N 82 READ (4,3) ( (IBND (L, J+K) ,L=1 ,2) , J=1, N) 83 K = K + N 84 200 CONTINUE 85 DO 300 1 = 1 , 200 86 CALL. FREAD (7, «2R: « ,YB (I) ,XB (I) ) 87 300 CONTINUE 88 RETURN 89 1 FORMAT (6012) 90 2 FORMAT (12) 91 3 FORMAT (40(I2,A1)) 92 END \ APPPENDIX I? - CATERPILLAR MODEL CODE SUBROUTINE ZERO 132 1 SUBROUTINE ZERO 2 C 3 C CALLED ONCE TO INITIAL THE G ARRAY TO 0.0 4 C 5 C COMMON /CM/ GOES HERE 6 JP1 = JHAX +1 7 DO 100 I = 1, IMAX 8 DO 100 J = 1, JP1 9 DO .100 K = 1, 15 10 G(K,J,I) = 0.0 11 100 CONTINOE 12 RETURN 13 END t APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE PLOTI 133 1 SUBROUTINE PLOTI 2 C 3 C PLOT INITIALIZATION 4 C 5 C XB AND TB ARE THE COORDINATES OF THE 6 C ROTATED OUTLINE 7 C 8 C COMMON /CM/ GOES HERE 9 SCALE = 60.0/HAXPOP 10 CALL ROTATI(FLOAT(ITURN),90.0-ITILT) 11 DX = 4.0/(MAX-MIX) 12 DY .= 4.0/(MAY-MI¥) 13 IF (IABS (IOPLOT) .NE. .1) RETURN 14 DO 20 1 = 1 , 200 15 X = XB(I) 16 Y = YB(I) 17 Z = 0 18 CALL ROTAT(Y,X,Z) 19 XB (I) = (X-HIX) *DX 20 YB(I) = (Y-MIY) *DY 21 IF (IOPLOT .EQ.,-1) CALL PLOTB (XB (I) , YB (I) , 2) 22 20 CONTINUE 23 RETURN 24 1 . FORMAT (10X,F3.0,10X,F3.0) 25 END APPPENDIX IV - CATERPILLAR MODEL CODE SOBROOTINE PLOTB 134 1 SUBROUTINE PLOTB(XrY,IP) 2 C 3 C PLOT BRANCH ROUTINE 4 C CALLS ADAGE PLOT IF IOPLOT [ 0 5 C ELSE USES CALCOMP PLOTS 6 C 7 C COMMON /CM/ GOES HERE 8 IF (IOPLOT .GT..0) GO TO 20 9 CALL PLOTA (X,Y,IP) 10 RETURN 11 20 CONTINOE 12 CALL PLOT (X,Y,IP) 13 RETORN 14 END APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE PLOTA 135 1 SUBROUTINE PLOTA (X,Y,IPEN) 2 C 3 C ADAGE PLOT ROUTINE 4 C 5 C COMMON /CM/ GOES HERE 6 LOGICAL*1 PENS(6000) 7 LOGICAL FIRST /.TRUE./ 8 DIMENSION XS(6000), YS (6000) 9 DATA MAXVEC /6000/ 10 DATA ISPOT,NSPOT / 2*1 / 11 IF( .NOT. FIRST ) GO TO 35 12 IF (NSPOT .NE. 1) GO TO 25 13 CALL AGTCON (XS(1),4,5,0) 14 CALL AGTCON (XS (2),6,3,0) 15 CALL AGTCON (XS (3) , 8,6, 0) 16 CALL AGTCON (XS (4) , 10, 2, 0) 17 XS (5) =0 18 NSPOT = 6 19 CALL AGTCON (LITON,34,1,0) 20 CALL AGTCON (LITOFF,36,1,0) 21 IPAGT = 0 22 25 CONTINOE 23 IF (IOPLOT .EQ.,-1) GO TO 30 24 CALL ADAGEB (XS,NSPOT-1,1,.TRUE.,2) 25 FIRST = .FALSE. 26 GO TO 35 27 30 CONTINUE 28 IF (IPEN .EQ. 3) GO TO 32 29 CALL AGTCVT (XS (NSPOT),X*2.0,Y*2.0,IPAGT,0) 30 NSPOT = NSPOT +1 31 IPAGT = 1 32 RETURN 33 32 CONTINUE 34 CALL ADAGEB (XS,NSPOT-1,1,.TRUE.,2) 35 FIRST = .FALSE. 36 35 CONTINUE 37 C*****CHECK FOR A PEN OF -3 ... ALL THROUGH ***** 38 IF( IPEN .LE. -3 ) GO TO 400 39 IF( ISPOT .LT. MAXVEC ) GO TO 200 40 PRINT 62 , IDAY, IYEAR 41 62 FORMAT ( »0 VECTOR OVERFLOS IDAY, IYEAR = «, 13, 10X, 13) 42 STOP 43 200 CONTINOE 44 XS ( ISPOT ) = X 45 YS ( ISPOT ) •=. Y 46 PENS( ISPOT ) = IPEN .EQ. 2 47 ISPOT = ISPOT +1 48 RETURN 49 400 ISPOT = ISPOT - 1 50 IF( ISPOT .LE. 0 ) RETURN 51 DO 450 1 = 1 , ISPOT 52 C******SET THE PEN BIT ************* APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE PLOTA 136 53 IPAGT = 0 54 IF( PENS( I ) ) IPAGT = 1 55 XAGT = XS(I)*2.0 56 CALL AGTCVT( XS (I) , XAGT, YS(I)*2.0, IPAGT, 0) 57 450 CONTINUE 58 IF(ISPOT .LT. 6000 ) ISPOT=ISPOT+1 59 CALL AGTCVT( XS(ISPOT), 0, 0 # 0, 1 ) 60 CALL ADAGEB ( XS , ISPOT , NSPOT , .FALSE., 2) 61 . CALL.IAITB(100) 62 CALL ADAGEB (LITON,1,5,.FALSE.,2) 63 CALL 8AITB(150) 64 CALL ADAGEB (LITOFF,1,5,.FALSE.,2) 65 CALL 8AITB(50) 66 ISPOT = 1 67 RETURN 68 END APPPENDIX I¥ - CATERPILLAB MODEL CODE SUBROUTINE ADAGEB 137 1 SUBROUTINE ADAGEB (I,J,K,L,M) 2 C 3 C ADAGE BRANCH ROUTINE 4 C SAVES THE DATA INSTEAD OF DISPLAYING IT 5 C IF IOSAVE NON-ZERO 6 C 7 C COMMON /CM/ GOES HERE 8 DIMENSION I (J) 9 IF (IOSAVE .NE. .0) GO TO 100 10 CALL AGTDSP <I,J,K,L,H) 11 RETURN 12 100 CONTINUE 13 WRITE (8) M,L#K,JrI 14 RETURN 15 END APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE WAITB 138 1 SUBROUTINE WAITB (H) 2 C 3 C WAITING BRANCH ROUTINE 4 C CALLS FOR A WAIT DURING DISPLAY UNLESS 5 C THE DATA IS BEING SAVED FOR DISPLAY 6 C LATER 7 C 8 C COMMON /CM/ GOES HERE 9 DATA I,J,K,L,H /1,1,-1,1,1/ 10 IF (IOSAVE .NE.,0) GO TO 100 11 CALL HTWAIT (N) 12 RETURN 13 100 CONTINUE 14 K = —N 15 WRITE (8) I,J,K,L,M 16 RETURN 17 END APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE INITAL 139 1 SUBROUTINE INITAL 2 C 3 C INITIALIZES THE G ARRAY 4 C 5 C IHEAD AND JHEAD POINT TO THE FIRST 6 C NON-ZERO GRID SQUARE IN G. THEN EACH 7 C NON-ZERO GRID SQUARE POINTS TO THE 8 C NEXT NON-ZERO ONE AND THE PREVIOUS ONE. 9 C THIS CREATES A DOUBLE LINKED LIST 10 C OF NON-ZERO GRID SQUARES. 11 C 12 C COMMON /CM/ GOES HERE 13 IRST = 0 14 IHEAD = 0 15 JHEAD =? 0 16 DO 40 I = 1, IMAX 17 DO 20 J = 1, JMAX 18 IF (G(5,J,I) .EQ. ,0) GO TO 20 19 IF (IRST .NE. 0) GO TO 10 20 G (12,J,I) = 0.0 21 G (13, J,I) = 0.0 22 IS = I 23 JS = J 24 IHEAD = IS 25 JHEAD = JS 26 IRST = 1 27 GO TO 20 28 10 CONTINOE 29 IG(14,JS,IS) = I 30 IG (15,JS,IS) = J 31 IG(12,J,I) = IS 32 IG (13,4*1) = JS 33 IS = I 34 JS = J 35 20 CONTINOE 36 40 CONTINOE 37 IG (14,JS,IS) =0 38 IG (15,JS,IS) -= 0 39 IF (IYEAR .NE. ,1) IDAY = 1 40 IF (IHEAD .EQ. 0) CALL HANGUP 41 RETURN 42 END APPPENDIX IV - CATERPILLAR MODEL CODE 1 SUBROUTINE HANGUP SUBROUTINE HANGUP* C C CALLED IF THE CATERPILLARS BECOME EXTINCT C C COMMON /CM/ GOES HERE WRITE (6,1) IDAY,IYEAR STOP 1 FORMAT (30H DISASTER...THEY ALL DIED ON I2,14HTH DAY IN THE , 1 . 12,10HTH YEAR!!!) END APPPENDIX IV - CATERPILLAR MODEL CODE SOBROOTINE OOTPOT 141 1 SOBROOTINE OOTPOT 2 C 3 C DEPENDS ON WHICH TYPE OF OOTPOT IS DESIRED 4 C 5 C COHMON /CM/ GOES HERE 6 LOGICAL L1, L2 7 EQUIVALENCE (L1,I0TYPE), (L2,MASK) 8 DO 100 II = 1, 5 9 MASK = 2**(II-1) 10 L2 = L1 .AND. .L2 11 IF (MASK .EQ..0) GO TO 100 12 CALL TOTALS (II) 13 IF (IOPLOT .NE. 0) GO TO 40 14 CALL PRNTER (II) 15 GO TO 100 16 40 CONTINOE 17 CALL PLOTER (II) 18 100 CONTINOE 19 RETURN 20 END APPPENDIX IV - CATERPILLAR MODEL CODE SOBROOTINE TOTALS 142 1 SOBROOTINE TOTALS (II) 2 C 3 C OOTPOTS POPULATION TOTALS ON THE PRINTER 4 C 5 C COMMON /CM/ GOES HERE 6 PMAX = 0 7 GTOT = 0 8 I = IHEAD 9 J = JHEAD 10 20 CONTINOE 11 GTOT = GTOT + G(II,J,I) 12 IF (G(II,J,I) . GT. , PMAX) PMAX = G(II,J,I) 13 IS = IG(14,J,I) 14 IF (IS .EQ. ,0) GO TO 25 15 J = IG(15,J,I) 16 I = IS 17 GO TO 20 18 25 CONTINOE 19 IF (IOTOTL .EQ.,0) GO TO 100 20 IF (IOPLOT .EQ.,0) GO TO 30 21 WRITE (6,2) IYEAR,IDAY, PMAX,GTOT 22 GO TO 100 23 30 CONTINOE 24 WRITE (6,1) HEAR, IDAY, PMAX, GTOT 25 100 CONTINUE 26 RETURN 27 1 FORMAT (30H1STATE OF SYSTEM IN YEAR 12X,I2,10H DAY ,12/ 28 1 30H POP MAX USED FOR SCALING E14.7/ 29 2 30H TOTAL POPOLATION E14.7) 30 2 FORMAT (30H-STATE OF SYSTEM IN YEAR 12X,I2,10H DAY ,12/ 31 1 30H POP MAX OSED FOR SCALING E14.7/ 32 2 30H TOTAL POPULATION E14.7) 33 END APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE PRNTER 143 1 SUBROUTINE PRNTER (II) 2 C 3 C OUTPUTS POPULATION DISTRIBUTIONS ON THE 4 C PRINTER 5 C 6 C COMMON /CM/ GOES HERE 7 INTEGER*2 ISYM, LEN, LINE 8 DIMENSION LINE (63), ISYM (36) 9 DATA LINE /63*1H / 10 DATA LEN /126/ 11 DATA ISYM /1H ,1H1,1H2,1H3,1H4,1H5,1H6,1H7, 1H8,1H9, 12 1 1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH, 1HI,1HJ, 13 2 1HK,1HL,1HH,1HN,1H0,1HP,1HQ,1HB, 1HS,1HT, 14 3 1HD,1HV,1HW,1HX,1HY,1HZ/ 15 IF (PMAX .LE. 0) GO TO 100 16 PMAX = 35/PMAX 17 M = 0 18 DO 40 I = 1, IMAX 19 DO 30 J .= 1, JHAX 20 IS = G(II,J,I)*PMAX • 1.99999 21 LINE(J+2) = ISYM (IS) 22 30 CONTINUE 23 LINE (2) = LINE(1) 24 KHAX = NBD(I) 25 DO 35 K = 1, KHAX 26 M = M + 1 27 L = IBND(1,M) • 1 28 IF (LINE (L) .NE.,LINE(1) .AND. „ L .LE.., JHAX) WRITE (6,1) 29 LINE (L) = IBND(2,M) 30 35 CONTINUE 31 CALL WRITE (LINE,LEN,0,LN0H,6) 32 40 CONTINUE 33 100 CONTINUE 34 RETURN 35 1 FORMAT (20H BOUNDRY ERROR ) 36 RETURN 37 END APPPENDIX IV - CATERPILLAR MODEL CODE SOBROOTINE PLOTER 144 SOBROOTINE PLOTER (II) C C C C C C C . C C C C . C C C C OOTPOTS POPOLATION DISTRIBUTIONS ON THE PLOTTER IHIST = 0 IHIST = 1 IHIST = 2 JHIST = 0 JHIST = 1 JHIST = 2 XO XO PREVIOOS POINT =0 PREVIOUS POINT =0 TWO .PREVIOUS POINTS PREVIOUS POINT ON GRID PREVIOUS POINT OFF GRID TWO PREVIOUS POINTS OFF GRID PREVIOUS X COORDINATE PREVIOUS Y COORDINATE =0 C COHMON /CH/ GOES HERE CALL HSPLOT (0.0/0.0,0,CLEAR) IF (IOPLOT .NE. ,1) GO TO 3 IIP = 3 DO 2 I = 1, 200 CALL PLOTB (XB(I) ,YB(I) ,IIP) IIP = 2 2 CONTINOE 3 CONTINUE K = 1 DO 14 I = 1, IMAX IHIST = 0 JHIST = 0 L = NBD(I) • K - 1 JS= IBND(1,K) - 1 JF= IBND(1,L) - 1 IP = -1 DO 13 J = JS, JF IF (J .NE.,0) GO TO 4 Z = 0 IHIST =1 JHIST = 1 GO TO 125 4 CONTINOE IF (IHIST .NE.,2) GO TO 6 IF (IABS(IOPLOT) .EQ. 1) GO TO 5 IF (G(10,J,I) ..NE. .0) GO TO 5 IHIST = 0 JHIST = 1 GO TO 12 5 CONTINOE IF (G (II, J,I) .NE. 0) GO TO 55 XO = J YO = I GO TO 13 55 CONTINOE CALL ROTAT(YO,XO,ZO) APPPENDIX IV - CATERPILLAR HODEL CODE SUBROUTINE PLOTER 145 54 X .= (XO-HIX) *DX 55 Y = (YO-HIY) *DY 56 IIP = -2 57 IF (IABS (IOPLOT) .EQ. 1) IIP = -1 58 CALL HSPLOT(X,Y,IIP,CLEAR) 59 IHIST = 0 60 GO TO 12 61 6 CONTINUE 62 IF (IABS (IOPLOT) .EQ. 1) GO TO 10 63 IF (JHIST .NE. 2) GO TO 8 64 IF (G(10,J,I) .NE.O) GO TO 7 65 XO = J 66 YO = I 67 ZO = 0 68 GO TO 13 69 7 CONTINUE 70 CALL ROTAT(YO,XO,ZO) 71 X = (XO-MIX) *DX 72 Y = (YO-HIY) *DY 73 IIP = -1 74 CALL ,HSPLOT(X,Y,IIP,CLEAR) 75 JHIST = 0 76 IHIST = 1 77 GO TO 10 78 8 CONTINUE 79 IF (G(10,J,I) .EQ. 0) GO TO 85 80 JHIST = 0 81 GO TO 10 82 85 CONTINUE 83 IF (JHIST .EQ-,0) GO TO 9 84 XO = J 85 YO = I 86 ZO = 0 87 JHIST = 2 88 GO TO 13 89 9 CONTINUE 90 JHIST = 1 91 GO TO 12 92 10 CONTINUE 93 IF (G(II,J,I) .EQ. ,0) GO TO 105 94 IHIST = 0 95 GO TO 12 96 105 CONTINUE 97 IF (IHIST .EQ.,0) GO TO 11 98 XO = J 99 YO = I 100 ZO = 0 101 IHIST = 2 102 GO TO 13 103 11 CONTINUE 104 IHIST = 1 105 12 CONTINUE 106 Z = 6(11,J,I) 107 IF (Z .NE. 0) Z = ALOG10 (Z) *SCALE APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE PLOTER 146 108 125 CONTINUE 109 X = J 110 Y = I 111 CALL ROTAT(Y,X#Z) 112 X = (X-MIX)*DX 113 Y = (Y-MIY)*D¥ 114 CALL HSPLOT (X,Y,IP#CLEAR) 115 IP = -2 116 13 CONTINUE 117 K = L + 1 118 14 CONTINUE 119 CALL PLOTB(12.0,O.0,-3) 120 RETURN 121 END APPPENDIX 17 - CATERPILLAR MODEL CODE SUBROUTINE ROTAT 147 1 SUBROUTINE ROTAT( X , Y , Z ) 2 C 3 C ROTATES 3-D VECTOR TO PROPER PERSPECTIVE 4 C 5 c**************************** 6 C* R O T AT * 7 C* 8 C* ROTATE THE X,Y,Z COORDINATES 9 C* * 11 XT = X 12 YT = ¥ 13 ZT = Z 14 X = XT*XC1 - YT*XC2 +ZT*XC3 15 Y = XT*YC1 + YT*YC2 - ZT*YC3 16 Z = XT*ZC1 • YT*ZC2 + ZT*ZC3 17 RETORN 18 ENTRY ROTATI( TORN, TILT ) 19 555 ASPECT = 0.0 20 DATA RC / .1745329E-1 / 21. T = TORN * RC 22 STUHN = SIN( T ) 23 CTURN = COS( T ) 24 T = TILT * RC 25 STILT = SIN( T ) 26 CTILT = COS( T ) 27 T = ASPECT * RC 28 SASP = SIN( T ) 29 CASP = COS{ T ) 30 C**** CALCULATE CONSTANTS 31 XC1 CTDRN*CTILT 32 XC2 = STURN*CTILT 33 XC3 = STILT 34 YC1 STURN*CASP + CTURN*STILT*SASP 35 YC2 CTURN*CASP - STURN*STILT*SASP 36 YC3 CTILT*SASP 37 ZC1 = STURN*SASP - CTURN*STILT*CASP 38 ZC2 = CTURN*SASP • STURN*STILT*CASP 39 ZC3 = CTILT*CASP 40 RETORN 41 END APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE HSPLOT 148 SUBROUTINE HSPLOT( XX , YY: , IPENN, CLEAR ) c * ^ * * * ***************************************** C* SUBROUTINE TO FACILITATE THE PLOTTING OF SURFACES * C* DIAGRAMS. HSPLOT HILL NOT PERMIT A LINE TO BE DRAHN AT A YHEIGHT* C* LESS THAN THE YHEIGHT REGISTERED IN THE YHT ARRAY FOR ANY X * C* POSITION. * C* * C* THUS IF A SURFACE IS DRAHN HITH THE AREA NEAREST THE * C*OBS£RVER BEING DRAHN FIRST HIDDEN LINES HILL BE REMOVED * C*FROH THE SURFACE BY THIS ROUTINE * C* * C* THE OUTPUT IS TO ROUTINE »PLOT» AND IS DESIGNED TO INTERFACE * C*DIRECTLY TO THE PLOT ROUTINES. ON OUTPUT: C* IPEN = 2 IS PEN DOHN * C* = 3 IS PEN DP. C* * C* RESTRICTION: PLEASE - NO PLOT LONGER THEN 30 INCHES. * C* * C* IF.,, * C* IPEN = +1 . +2 , +3 THE MOVEMENT IS MADE HITHOOT CHECKING. * C* IPEN = 0 THE YHT ARRAY IS CLEARED TO -0.5. * C* IPEN = -1 THE PEN IS MOVED VIRTUALLY. IF THE POSITION INDICATED * C* IS IN THE CLEAR THE PEN IS MOVED THERE OP. IF THE * C* POSITION IS HIDDEN THE PEN IS NOT ACTOALLY MOVED * C* BUT THE POSITON INDICATORS ARE UPDATED., ' * C* IPEN = -2 THE MOVEMENT IS MADE HITH CHECKING AND THE YHT ARRAY C* IS UPDATED.. C* IPEN = -3 THE MOVEMENT IS HADE HITH APPPENDIX IV - CATERPILLAR HODEL CODE SUBROUTINE HSPLOT CHECKING BUT THE YHT ARRAY C* IS NOT UPDATED * C* IPEN = -H THE MOVEMENT IS ONLY SIMULATED., THE YHT ARRAY IS * C* UPDATED BUT THE PEN IS NOT ACTUALLY HOVED. , C* * C* CLEAR - HILL BE RETURNED .TRUE. IF THE PEN HAS ENTIRELY IN THE * C* CLEAR (NOT HIDDEN IN ANY HAY) . IT HILL BE .FALSE. ,OTHERHISE* C* IT IS USED PRIMARILY TO OPTIMIZE PEN UP AND DOHN MOTION BY * C* ELIMINATING THE NEED TO UNNECCESARILY LIFT THE PEN IF A * C* LINE IS TO BE RETRACED. * C* C* * DIMENSION YHT( 3201 ) C C IF PEN = .TRUE. THE PEN IS IN AN AREA WHERE PLOTTING MAY OCCUR C LOGICAL CLEAR LOGICAL PEN , FLAG4 , FLAG3 DATA NERROR/ 0 / X = XX ¥ = YY \ IPEN = IPENN IF( IPEN .NE. 0 ) GO TO 100 DO 50 I = 1 , 3201 50 YHT ( I ) = ( -0.5 ) PEN =- .TRUE., RETURN 100 IF( IPEN .LT. 0 ) GO TO 200 CALL PLOTB( X , Y , IPEN ) XO = X YO = Y IX = X * 100.0 • 0.5 IX = IX + 200 PEN = .TRUE. IF( YHT( IX ) .GT. Y • 0.005 ) PEN = .FALSE. RETURN 200 CONTINUE IF (IPEN .NE.,( -1 ) ) GO TO 201 XO = X YO = Y APPPENDIX IV - CATERPILLAR HODEL CODE SUBROUTINE HSPLOT 150 69 IX = X * 100.0 • 0.5 70 IX = IX + 200 71 PEN .= • .FALSE. , 72 IF( YHT( IX ) .GT. Y • 0.005) RETORN 73 PEN = .TROE. 74 CALL PLOTB( X , Y , 3 ) 75 RETORN 76 201 . CONTINUE 77 C****** 'CLEAR* SILL BE TRUE IF THE LINE IS NOT HIDDEN ON ANY ******** 78 c****** PORTION OF THE JOURNEY **************************** 79 CLEAR = PEN 80 C***** RANGE OF THE X HOVEHENT 81 IF( X .LT.„ ( -1.9 ) .OR. X .GT. ,30.0 ) GO TO 99 82 IX = X * 100.0 • 0.5 83 IXO = XO * 100.0 • 0.5 84 DX = IABS( IX - IXO ) 85 C**** THE Y INCREMENT 86 DYI = I..- YO 87 IF( DX .NE. 0.0 ) DYI = DYI / DX 88 ADYI = ABS ( DYI ) 89 IXEND = IX • 200 90 IXLOC = IXO 200 91 INC =1 92 IF( IXLOC .GT. IXEND ) INC = ( -1 ) 93 IXE = IXEND + INC 94 YS = YO 95 C 96 C PREAMBLE -97 C CHECK FOR NO X MOVEMENT. , 98 C 99 IF( DX .NE. ,0.0 ) GO TO 340 100 YS = Y 101 YSOLD = YS 102 IXLO = IXLOC 103 C 104 C************* S T A R T OF M A I N L O O p ***************** 105 C 106 300 CONTINUE 107 IF( YHT( IXLOC ).GT. YS • 0.005 ) GO TO 330 108 C 109 C PEN IS IN THE CLEAR SO UPDATE THE . •YHT* ARRAY AND MOVE 110 C THE PEN UP TO THE LAST POSITION HITH THE PEN UP IF THE LAST 111 C POSITION HAS HIDDEN. 112 C 113 IF( PEN ) GO TO 350 114 XLOC = FLOAT( IXLO - 200 ) * 0.01 115 C 116 C USE THE Y-VALUE THAT HILL CAUSE APPPENDIX I? - CATERPILLAR MODEL CODE SUBROUTINE HSPLOT 151 THE SHALLEST VARIATION 117 C IN THE X DIRECTION. 118 C 119 YPLOT = YSOLD 120 IF( ADYI.GT. ABS(YHT(IXLO) - YHT( IXLOC) ) ) YPLOT = YHT (IXLO ) 121 CALL PLOTB( XLOC , YPLOT , +3 ) 122 PEN = .TRUE. 123 350 YHT ( IXLOC ) = YS 124 GO TO 340 125 330 CONTINUE 126 C 127 C THE PEN IS HIDDEN - CAN*T PLOT IN THIS AREA. IF THE 128 C LAST POSITION HAS IN THE CLEAR MOVE THE PEN THERE DOHN. 129 C 130 IF( .NOT.,PEN ) GO TO 340 131 XLOC = FLOAT( IXLO - 200 ) * 0.01 132 YPLOT = YS 133 IF( ADYI.GT. ABS(YHT(IXLO) - YHT( IXLOC) ) ) YPLOT = YHT (IXLOC) 134 CALL PLOTB( XLOC , YPLOT , +2 ) 135 PEN = .FALSE. 136 CLEAR = .FALSE. 137 340 CONTINUE 138 YSOLD = YS 139 YS = YS + DYI 140 IXLO = IXLOC 141 IXLOC = IXLOC • INC 142 IF( IXLOC .NE. IXE ) GO TO 300 143 C 144 C************** E N D OF M A I N L O O P ********************* 145 C 146 C SAVE THE LAST POSITION. DON'T ACTUALLY MOVE THE PEN 147 C IF IT IS HIDDEN.. 148 C 149 XO = X 150 YO = Y 151 IF( .NOT. PEN ) RETURN 152 CALL PLOTB( X , Y , 2 ) 153 RETURN 154 , 99 CONTINUE 155 IF( HERROR .GT..0 ) RETURN 156 NERROR = NERROR + 1 . 157 HRITE( 6 , 991 ) X 158 991 FORMAT( 1H0 , 51HERROR IN CALL TO •HSPLOT*. ATTEMPT TO PLOT AT 159 1= , G16.8 , 9H INCHES. / ) 160 RETURN 161 END APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE STEP 152 1 SUBROUTINE STEP 2 C 3 C GOES THROUGH LINK LIST OF G«S ONCE TO 4 C INTEGRATE AHEAD ONE STEP OF NDAYS. 5 C 6 C IHPT IS THE POINTER TO THE HEATHER FOR THIS 7 C DAY 8 C 9 C COMMON /CM/ GOES HERE 10 10 CONTINUE 11 IF (IYEAR .LT., IHETH (2,IHPT) ) GO TO 40 12 IF (HEAR . EQ. IHETH (2 ,IHPT) ) GO TO 20 13 IHPT = IHPT + 1 14 GO TO 10 15 20 CONTINUE 16 IF (IDAY .LE.,IHETH(1,IHPT)) GO TO 40 17 IHPT = IHPT • 1 18 GO TO 20 19 40 CONTINUE 20 I = IHEAD 21 J = JHEAD 22 50 CONTINUE 23 CALL INTGRT (G (1, J,I) ,G (1, J,I) ,IWETH (3, IHPT) ) 24 IF (G(5,J,I) .NE. 0) GO TO 60 25 IF (IG(12,J,I) .NE. 0) GO TO 55 26 IHEAD = IG(14,J,I) ^ 27 JHEAD = IG(15,J,I) 28 IF (IHEAD .EQ.,0) CALL HANGUP 29 IG (12, JHEAD,IHEAD) = 0 30 IG (13,JHEAD,IHEAD) = 0 31 GO TO 60 32 55 CONTINUE 33 IP = IG(12,J,I) 34 JP = IG(13,J,I) 35 IG(14,JP,IP) = IG(14,J,I) 36 IG(15,JP,IP) = IG(15,J,I) 37 IP = IG(14,J,I) 38 IF (IP .EQ. 0) GO TO 80 39 JP = IG(15,J,I) 40 IG(12,JP,IP) = IG(12,J,I) 41 IG(13,JP,IP) = IG(13,J,I) 42 60 CONTINUE 43 1 1 = IG(14,J,I) 44 J = IG (15,J,I) 45 I = II 46 IF (I .NE. 0) GO TO 50 47 80 CONTINUE 48 RETURN 49 END APPPENDIX I? - CATERPILLAR HODEL CODE SUBROUTINE INTGRT 153 1 SUBROUTINE INTGRT (GS,IGS#I8) 2 C 3 C INTEGRATES ONE GRID SQUARE NDATS AHEAD IN 4 C TIHE 5 . C 6 ' C COHHON /CH/ GOES HERE 7 DIMENSION GS(12), IGS(12) 8 PCT1 = GS(1)/GS(5) - .2 9 IF (PCT1 .GT. 0) PCT1 = 0 10 LOCH = IH*IGS(10) 11 GS (5) = 0 12 DO 20 1 = 1 , 4 13 GS(I+5) = (C1(I) • C2(I)*IDAX + C3 (I) *LOCH • CH (I) *PCT1) *GS (I) 14 GS(I) = GS(I) + GS(I+5)*NDAY 15 IF (GS(I) .LT. 1) GS(I) = 0 16 GS(5) = GS(5) • GS(I) 17 20 CONTINUE 18 RETURN 19 END APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE DSPRS 154 1 SUBROUTINE DSPRS 2 C 3 C DISPERSES THE ADULT HOTHS TO NEW GRID 4 C SQUARES AND UPDATES FOR NE8 EGG MASSES 5 C 6 C COMMON /CM/ GOES HERE 7 DO 4 I = 1, IHAX 8 DO .4 J = 1, JMAX 9 DO 4 K = 5, 9 10 G (K,J,I) = 0 11 4 CONTINOE 12 I = IHEAD 13 J = JHEAD 14 C 15 C 16 C 17 IASUB = RAND (1.0) 18 5 CONTINOE 19 X = I - .5 20 Y = J - .5 21 DO 100 K = 1, 3 22 L = G(K,J,I) 23 IF (L .EQ- 0) GO TO 100 24 DO 95 H = 1, L 25 IASUB = FHAND(1.0) *108.0 + 1.0 26 III = I 27 JJJ = J 28 IF (IASOB .LT. 82) GO TO 10 29 INC = 1 30 JNC = -1 31 II = I 32 JJ = J - 1 33 N = 110 - IASOB 34 S -STAB (N) 35 C CTAB (N) 36 GO TO 40 37 10 CONTINOE 38 IF (IASOB .LT. 56) GO TO 20 39 INC =. -1 40 JNC = -1 41 II = I - 1 4 2 JJ = J - 1 43 . N = IASOB - 54 44 S -STAB (N) 45 C -CTAB (N) 46 GO TO 40 47 20 CONTINUE 48 IF (IASUB .LE. ,28) GO TO 30 49 INC = -1 50 JNC = 1 51 II = I - 1 52 JJ = J 53 N 56 - IASOB APPPENDIX IV - CATERPILLAR MODEL CODE SOBROOTINE DSPRS 155 54 S STAB (N) 55 C -CTAB (N) 56 GO TO 40 57 30 CONTINOE 58 INC = 1 59 JNC = 1 60 II = I 61 JJ = J 62 S STAB (IASOB) 63 C CTAB (IASOB) 64 40 CONTINOE 65 DT = FRAND(2.0)*DTAB(K) 66 D1 = (JJ-¥)/S 67 D2 = (II-X)/C 68 DOLD= 0 ) 69 D = 0 70 50 CONTINOE 71 DIF = D1 - D2 72 IF (ABS (DIF) .GT. 1.0E-3) GO TO 60 73 D = D • (D1-DOLD)*IABS(IG (11,JJJ,III)) 74 II = II + INC 75 JJ = JJ • JNC 76 III = III* INC 77 JJJ = JJJ* JNC 78 DOLD= D1 79 D1 = (JJ-IJ/S 80 D2 = (II-X)/C 81 GO TO 80 82 60 CONTINOE 83 IF (DIF .LT. ,0) GO TO 70 84 D = D + (D2-DOLD)*IABS(IG(11,JJJ,III) ) 85 II = II + INC 86 III = III* INC 87 DOLD= D2 88 D2 = (II-X)/C 89 GO TO 80 90 70 CONTINOE 91 D = D • (D1-DOLD)*IABS(IG(11,JJJ,III) ) 92 JJ = JJ • JNC 93 JJJ = JJJ+ JNC 94 DOLD= D1 95 D1 = (JJ-Y)/S 96 80 CONTINOE 97 IF (D .GT. DT) GO TO 85 98 IF (II .GT. IMAX) GO TO 95 99 IF (JJ .GT. JMAX) GO TO 95 100 IF (II .LT. 1). GO TO 95 101 IF (JJ .LT. 1) GO TO 95 102 GO TO 50 103 85 CONTINOE 104 IF (IG(11,JJJ,III) • EQ. -1) GO TO 95 105 DO 90 N = 5, 9 106 G(N,JJJ,III) = G(N,JJJ,III) • ETAB (N-4 ,K) 107 90 CONTINOE APPPENDIX IV - CATERPILLAR MODEL CODE SUBROUTINE DSPRS 156 108 95 CONTINUE 109 100 CONTINUE 110 II = IG(14,J,I) 111 J = IG(15,J,I) 112 I = II 113 IF (I .GT. 0) GO TO 5 114 C 1 115 C 116 C 117 DO 500 I = 1, IMAX 118 DO 500 J = 1, JMAX 119 DO 500 K .= 1, 4 120 G(K,J#I) = G(K+5,J#I) 121 500 CONTINUE 122 RETURN 123 END
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Of caterpillars and models
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Of caterpillars and models Steinhoff, Richard Terrell 1975
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Of caterpillars and models |
Creator |
Steinhoff, Richard Terrell |
Publisher | University of British Columbia |
Date Issued | 1975 |
Description | An approach to modeling in general, its background, and its implementation can be described by examining its application to a particular system. The Western Tent Caterpillar is a productive subject of study because there is extensive data available for determining the relationships within the system. Individual differences between caterpillars and the weather are the dominant factors involved. They are incorporated into a set of ordinary differential equations that express the rate of change in the caterpillar population over time. These equations are numerically integrated to calculate the population at some terminal time. The results over a six year period agree reasonably well with observed field data. The three dimensional plot of the population gives a visual appreciation for the distribution over time. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-28 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0051765 |
URI | http://hdl.handle.net/2429/19241 |
Degree |
Master of Science - MSc |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1975_A6_7 S74.pdf [ 6.57MB ]
- Metadata
- JSON: 831-1.0051765.json
- JSON-LD: 831-1.0051765-ld.json
- RDF/XML (Pretty): 831-1.0051765-rdf.xml
- RDF/JSON: 831-1.0051765-rdf.json
- Turtle: 831-1.0051765-turtle.txt
- N-Triples: 831-1.0051765-rdf-ntriples.txt
- Original Record: 831-1.0051765-source.json
- Full Text
- 831-1.0051765-fulltext.txt
- Citation
- 831-1.0051765.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0051765/manifest