f1-1

1.1ProgramOrganizationandControlStructures 5Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Previous Routines Omitted from ThisEdition Name(s) Replacement(s) Comment ADI mglin ormgfas better method COSFT cosft1 orcosft2 choice of boundary conditions CEL,EL2 rf ,rd,rj,rc better algorithms DES,DESKS ran4 now uses psdes was too slow MDIAN1,MDIAN2 select ,selip more general QCKSRT sort name change ( SORTisnow hpsort) RKQC rkqs better method SMOOFT useconvlvwith coefficients from savgol SPARSE linbcg more general is sometimes quite difficult. We have not attempted this, and we do not pretend to any degree of bibliographical completeness in this book. For topics where a substantial secondary literature exists (discussion in textbooks, reviews, etc.) we have consciously limited our references to a few of the more useful secondarysources, especially those with good references to the primary literature. Where the existing secondary literature is insufficient, we give references to a few primary sources that are intended to serve as starting points for further reading, not as complete bibliographies for the field. Theorderinwhichreferencesarelistedisnotnecessarilysignificant. Itreflectsa compromisebetweenlistingcitedreferencesintheordercited,andlistingsuggestions for furtherreadingin a roughlyprioritizedorder,with the most useful ones first. The remaining two sections of this chapter review some basic concepts of programming (control structures, etc.) and of numerical analysis (roundoff error, etc.). Thereafter, we plunge into the substantive material of the book. CITED REFERENCES AND FURTHER READING: Meeus, J. 1982, Astronomical Formulae for Calculators , 2nd ed., revised and enlarged (Rich- mond, VA: Willmann-Bell). [1] 1.1 Program Organization and Control Structures Wesometimesliketopointoutthecloseanalogiesbetweencomputerprograms, on the one hand, and written poetry or written musical scores, on the other. All three present themselves as visual media, symbols on a two-dimensional page or computerscreen. Yet, in all three cases, the visual, two-dimensional, frozen-in-time representation communicates (or is supposed to communicate) something rather 6 Chapter1. PreliminariesSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).different,namelya processthat unfoldsin time . A poemis meantto beread;music, played; a program,executed as a sequential series of computer instructions. Inallthreecases,thetargetofthecommunication,initsvisualform,isahuman being. The goal is to transfer to him/her, as efficiently as can be accomplished, the greatest degree of understanding, in advance, of how the process willunfold in time. In poetry, this human target is the reader. In music, it is the performer. In programming, it is the program user. Now, you may object that the target of communication of a program is not a human but a computer, that the program user is only an irrelevant intermediary, a lackey who feeds the machine. This is perhaps the case in the situation wherethe business executive pops a diskette into a desktop computer and feeds that computer a black-box program in binary executable form. The computer, in this case,doesn’tmuchcarewhetherthatprogramwas writtenwith“goodprogrammingpractice” or not. Weenvision,however,thatyou,thereadersofthisbook,areinquiteadifferent situation. You need, or want, to know not just whata program does, but also how it does it, so that you can tinker with it and modify it to your particular application. You need others to be able to see what you have done, so that they can criticize or admire. In such cases, where the desired goal is maintainable orreusablecode, the targets of a program’s communication are surely human, not machine. One key to achieving good programming practice is to recognize that pro- gramming, music, and poetry — all three being symbolic constructs of the human brain — are naturally structured into hierarchies that have many different nested levels. Sounds(phonemes)formsmall meaningfulunits(morphemes)whichin turnformwords; words groupinto phrases, whichgroupinto sentences; sentences make paragraphs, and these are organized into higher levels of meaning. Notes form musical phrases, which form themes, counterpoints, harmonies, etc.; which form movements, which form concertos, symphonies, and so on. The structure in programs is equally hierarchical. Appropriately, good pro- gramming practice brings different techniques to bear on the different levels [1-3]. At a low level is the asciicharacter set. Then, constants, identifiers, operands, operators. Then program statements, like a(j+1)=b+c/3.0 . Here, the best pro- gramming advice is simply be clear, or (correspondingly) don’t be too tricky .Y o u might momentarily be proud of yourself at writing the single line k=(2-j)*(1+3*j)/2 if you want to permute cyclically one of the values j=( 0 ,1,2)into respectively k=( 1 ,2,0). You will regret it later, however, when you try to understand that line. Better, and likely also faster, is k=j+1 if (k.eq.3) k=0 Many programmingstylists would even argue for the ploddingly literal if (j.eq.0) then k=1 else if (j.eq.1) then k=2 1.1ProgramOrganizationandControlStructures 7Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).else if (j.eq.2) then k=0 else pause ’never get here’ endif onthegroundsthatitisbothclearandadditionallysafeguardedfromwrongassump- tions about the possible values of j. Our preference among the implementations is for the middle one. In this simple example, we have in fact traversed several levels of hierarchy: Statements frequently come in “groups” or “blocks” which make sense only taken as a whole. The middle fragment above is one example. Another is swap=a(j) a(j)=b(j)b(j)=swap whichmakes immediatesense to anyprogrammeras the exchangeof two variables, while sum=0.0 ans=0.0n=1 is verylikelyto bean initializationof variablespriortosome iterativeprocess. This levelofhierarchyinaprogramisusuallyevidenttotheeye. Itisgoodprogramming practiceto putin commentsat this level, e.g.,“initialize”or “exchangevariables.” The next level is that of control structures . These are things like the if...then ...elseclauses in the example above, doloops, and so on. This level is sufficiently important, and relevant to the hierarchical level of the routines in this book, that we will come back to it just below. At still higher levels in the hierarchy, we have (in FORTRAN) subroutines, functions, and the whole “global” organization of the computational task to be done. In the musical analogy, we are now at the level of movements and complete works. At these levels, modularization andencapsulation become important programming concepts, the general idea being that program units should interactwithoneanotheronlythroughclearlydefinedandnarrowlycircumscribedinterfaces. Good modularization practice is an essential prerequisite to the success of large, complicated software projects, especially those employing the efforts of more thanoneprogrammer. Itisalsogoodpractice(ifnotquiteasessential)inthelessmassive programmingtasks that an individualscientist, or readerof this book,encounters. Somecomputerlanguages,suchasModula-2and C++,promotegoodmodular- ization with higher-level language constructs, absent in FORTRAN-77. In Modula-2, for example, subroutines, type definitions, and data structures can be encapsulatedinto “modules” that communicate through declared public interfaces and whose internal workings are hidden from the rest of the program [4]. In the C++language, the key concept is “class,” a user-definablegeneralizationof data type that providesfor data hiding, automatic initialization of data, memory management, dynamic typing, and operatoroverloading(i.e., the user-definableextensionof operators like +and*so as to be appropriate to operands in any particular class) [5]. Properly used in defining the data structures that are passed between program units, classes 8 Chapter1. PreliminariesSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).can clarify and circumscribe these units’ public interfaces, reducing the chances of programming error and also allowing a considerable degree of compile-time andrun-time error checking. Beyond modularization, though depending on it, lie the concepts of object- orientedprogramming . Hereaprogramminglanguage,suchas C++orTurboPascal 5.5 [6], allows a module’spublicinterfaceto acceptredefinitionsoftypesor actions, and these redefinitions become shared all the way down through the module’s hierarchy(so-called polymorphism ). Forexample,aroutinewrittentoinvertamatrix ofrealnumberscould—dynamically,atruntime—bemadeabletohandlecomplex numbers by overloading complex data types and corresponding definitions of thearithmeticoperations. Additionalconceptsof inheritance (theabilitytodefineadata type that “inherits” all the structure of another type, plus additional structure of its own), and object extensibility (the ability to add functionality to a module without access to its source code, e.g., at run time), also come into play. We have not attempted to modularize, or make objects out of, the routines in this book, for at least two reasons. First, the chosen language, FORTRAN-77, does not really make this possible. Second, we envision that you, the reader, might want to incorporate the algorithms in this book, a few at a time, into modules or objects withastructureofyourownchoosing. Theredoesnotexist,atpresent,astandardor accepted set of “classes” for scientific object-oriented computing. While we might have tried to invent such a set, doing so would have inevitably tied the algorithmiccontent of the book (which is its raison d’ ˆetre) to some rather specific, and perhaps haphazard, set of choices regarding class definitions. On the other hand, we are not unfriendly to the goals of modular and object- oriented programming. Within the limits of FORTRAN, we have therefore tried to structureourprogramstobe“objectfriendly,”principallyviathecleardelineationof interface vs. implementation( §1.0) and the explicit declarationof variables. Within our implementation sections, we have paid particular attention to the practices of structured programming , as we now discuss. ControlStructures An executing program unfolds in time, but not strictly in the linear order in which the statements are written. Programstatements that affect the order in which statements are executed, or that affect whether statements are executed, are calledcontrolstatements . Controlstatementsnevermakeusefulsensebythemselves. They makesenseonlyinthecontextofthegroupsorblocksofstatementsthattheyinturn control. If you think of those blocks as paragraphs containing sentences, then the control statements are perhaps best thought of as the indentation of the paragraph and the punctuationbetweenthe sentences, not the words within the sentences. We can now say what the goal of structured programming is. It is to make programcontrol manifestly apparentin the visual presentationof the program .Y o u see that this goal has nothing at all to do with how the computer sees the program.Asalreadyremarked,computersdon’tcarewhetheryouusestructuredprogramming or not. Human readers, however, docare. You yourself will also care, once you discover how much easier it is to perfect and debug a well-structured program than one whose control structure is obscure. 1.1ProgramOrganizationandControlStructures 9Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).You accomplish the goals of structured programming in two complementary ways. First, you acquaint yourself with the small number of essential controlstructures that occur over and over again in programming, and that are therefore givenconvenientrepresentationsinmostprogramminglanguages. Youshouldlearn to think about your programmingtasks, insofar as possible, exclusively in terms ofthese standardcontrolstructures. In writing programs,youshouldget into the habit of representingthese standardcontrolstructuresin consistent,conventionalways. “Doesn’t this inhibit creativity?” our students sometimes ask. Yes, just as Mozart’s creativity was inhibited by the sonata form, or Shakespeare’s by the metrical requirementsof the sonnet. The pointis that creativity, whenit is meant tocommunicate,does wellundertheinhibitionsofappropriaterestrictionsonformat. Second, you avoid, insofar as possible, control statements whose controlled blocks or objects are difficult to discern at a glance. This means, in practice, thatyou must try to avoid statement labels and goto’s.It is not the goto’s that are dangerous (although they do interrupt one’s reading of a program); the statement labels are the hazard. In fact, whenever you encounter a statement label whilereading a program, you will soon become conditioned to get a sinking feeling in the pit of your stomach. Why? Because the following questions will, by habit, immediatelyspringtomind: Where didcontrolcome fromina branchtothis label? It could be anywhere in the routine! What circumstances resulted in a branch to this label? They could be anything! Certainty becomes uncertainty, understandingdissolves into a morass of possibilities. Some older languages, notably1966 FORTRAN and to a lesser extent FORTRAN- 77,requirestatementlabelsintheconstructionofcertainstandardcontrolstructures. We will see this in more detail below. This is a demerit for these languages. In such cases, you must use labels as required. But you should never branch to them independently of the standard control structure. If you must branch, let it be to an additionallabel,onethat is notmasqueradingas partofa standardcontrolstructure. We call labels that are part of a standard construction and never otherwise branchedto tame labels . They do not interfere with structured programmingin any way, except possibly typographically as distractions to the eye. Some examples are now in order to make these considerations more concrete (see Figure 1.1.1). Catalog of Standard Structures Iteration. InFORTRAN, simple iteration is performed with a doloop, for example do 10 j=2,1000 b(j)=a(j-1)a(j-1)=j 10 continue Notice how we always indent the block of code that is acted upon by the control structure, leaving the structure itself unindented. The statement label 10in this example is a tame label. The majority of modern implementations of FORTRAN-77 provide a nonstandard language extension that obviates the tame label. Originally 10 Chapter1. PreliminariesSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).yes no DO iteration (a)false true DO WHILE iteration (b) true false BREAK iteration (d)false true DO UNTIL iteration (c)iteration complete? block increment indexwhile condition until conditionblock break conditionblock blockblock Figure 1.1.1. Standard control structures used in structured programming: (a) DO iteration; (b) DO WHILE iteration; (c) DO UNTIL iteration; (d) BREAK iteration; (e) IF structure; (f) obsolete form ofDO iteration found in FORTRAN-66 , where the block is executed once even if the iteration condition is initially not satis fied. 1.1ProgramOrganizationandControlStructures 11Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).if condition blocktrueelse if condition blockfalse true. . . . . .false else blockelse if condition blockfalse true IF structure (e) iteration complete?increment index noblock FORTRAN-66 DO (obsolete) (f)yes Figure1.1.1. Standardcontrolstructuresusedinstructuredprogramming(seecaptiononpreviouspage). 12 Chapter1. PreliminariesSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).introduced in Digital Equipment Corporations ’s VAX-11 FORTRAN, the“enddo” statement is used as do j=2,1000 b(j)=a(j-1)a(j-1)=j enddo Infact,itwasaterriblemistakethattheAmericanNationalStandardfor FORTRAN-77 (ANSI X3.9 –1978) failed to provide an enddoor equivalent construction. This mistake by the people who write standards, whoever they are, presents us now,more than 15 years later, with a painful quandary: Do we stick to the standard, and clutterourprogramswith tamelabels? Ordowe adopta nonstandard(albeitwidely implemented) FORTRAN construction like enddo? We have adopteda compromiseposition. Standards,evenimperfectstandards, areterriblyimportantandhighlynecessaryinatimeofrapidevolutionincomputers andtheirapplications. Therefore,all machine-readableformsofourprograms(e.g., the diskettes that you can order from the publisher —see back of this book) are strictly FORTRAN-77 compliant. (Well, almoststrictly: there is a minor anomaly regardingbitmanipulationfunctions,seebelow.) Inparticular, doblocksalwaysend with labeled continue statements, as in the first example above. In the printed version of this book, however, we make use of typography to mitigatethestandard ’sdeficiencies. Thestatementlabelthatfollowsthe doisprinted in small type —as a signal that it is a tame label that you can safely ignore. And, theword“continue”is printedas “enddo”, whichyoumayregardas a verypeculiar change of font! The exampleabove, in our adoptedtypographicalformat, is do10j=2,1000 b(j)=a(j-1) a(j-1)=j enddo 10 (Noticethatwe alsotakethetypographicallibertyofwritingthetamelabelafterthe “continue”statement, rather than before.) A nested doloop looks like this: do12j=1,20 s(j)=0. do11k=5,10 s(j)=s(j)+a(j,k) enddo 11 enddo 12 Generally, the numerical values of the tame labels are chosen to put the enddo’s (labeled continue ’s onthediskette)intoascendingnumericalorder,hencethe do12 before the do11in the above example. IF structure. In this structure the FORTRAN-77 standard is exemplary. Here is a working program that consists dominantly of ifcontrol statements: 1.1ProgramOrganizationandControlStructures 13Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).FUNCTION julday(mm,id,iyyy) INTEGER julday,id,iyyy,mm,IGREG PARAMETER (IGREG=15+31*(10+12*1582)) Gregorian Calendar adopted Oct. 15, 1582. In this routine julday returns the Julian Day Number that begins at noon ofthe calendar date specified by month mm,d a y id,a n dy e a r iyyy, all integer variables. Positive year signifies A.D.; negative, B.C. Remember that the year after 1 B.C. was 1 A.D. INTEGER ja,jm,jyjy=iyyyif (jy.eq.0) pause ’julday: there is no year zero’ if (jy.lt.0) jy=jy+1 if (mm.gt.2) then Here is an example ofa block IF-structure. jm=mm+1 else jy=jy-1 jm=mm+13 endif julday=365*jy+int(0.25d0*jy+2000.d0)+int(30.6001d0*jm)+id+1718995 if (id+31*(mm+12*iyyy).ge.IGREG) then Test whether to change to Gregorian Calen- dar. ja=int(0.01d0*jy) julday=julday+2-ja+int(0.25d0*ja) endif returnEND (Astronomers number each 24-hour period, starting and ending at noon, with a unique integer, the Julian Day Number [7]. Julian Day Zero was a very long time ago; a convenient reference point is that Julian Day 2440000 began at noon of May 23, 1968. If you know the Julian Day Number that begins at noon of a given calendar date, then the day of the week of that date is obtained by adding 1and taking the result modulo base 7; a zero answer corresponds to Sunday, 1 to Monday, ..., 6 to Saturday.) Do-While iteration. Most good languages, except FORTRAN, provide for structures like the following Cexample: while (n<1000) { n=2*n; j++; InCthis has the meaning j=j+1 . } In fact, many FORTRAN implementationshave the nonstandardextension do while (n.lt.1000) n=2*n j=j+1 enddo Withinthe FORTRAN-77standard,however,thestructurerequiresatamelabel: 17if (n.lt.1000) then n=2*n j=j+1 goto 17 endif 14 Chapter1. PreliminariesSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).ThereareotherwaysofconstructingaDo-Whilein FORTRAN,butwetrytouse the aboveformat consistently. You will quickly get used to a statement like 17ifas signaling this structure. Notice that the two final statements are not indented, since they are part of the control structure, not of the inside block. Do-Until iteration. InPascal, for example, this is rendered as REPEAT n:=n DIV 2; Pascal’s integer divide is DIV. k:=k+1; UNTIL (n=1); InFORTRAN we write 19continue n=n/2 k=k+1 if (n.ne.1) goto 19 Break. In this case, you have a loop that is repeated inde finitely until some condition tested somewhere in the middle of the loop (and possibly tested in more than one place) becomes true. At that point you wish to exit the loop and proceedwith whatcomes afterit. Standard FORTRAN does notmakethis structureaccessible without labels. We will try to avoid using the structure when we can. Sometimes, however, it is plainly necessary. We do not have the patience to argue with thedesigners of computer languages over this point. In FORTRAN we write 13continue [statements before the test ] if ( ···) goto 14 [statements after the test ] goto 13 14continue Hereis aprogramthatusesseveraldifferentiterationstructures. Oneofuswas once asked, for a scavenger hunt, to find the date of a Friday the 13th on which the moon was full. This is a programwhich accomplishes that task, giving incidentallyall other Fridays the 13th as a by-product. PROGRAM badluk INTEGER ic,icon,idwk,ifrac,im,iybeg,iyend,iyyy,jd,jday,n, * julday REAL TIMZON,fracPARAMETER (TIMZON=-5./24.) Time zone −5is Eastern Standard Time. DATA iybeg,iyend /1900,2000/ The range ofdates to be searched. C USES flmoon,julday write (*,’(1x,a,i5,a,i5)’) ’Full moons on Friday the 13th from’, * iybeg,’ to’,iyend do12iyyy=iybeg,iyend Loop over each year, do11im=1,12 and each month. jday=julday(im,13,iyyy) Is the 13th a Friday? idwk=mod(jday+1,7) if(idwk.eq.5) then n=12.37*(iyyy-1900+(im-0.5)/12.) This value nis a first approximation to how many full moons have occurred since 1900. We will feed it into the phase routine and adjust it up or down until 1.1ProgramOrganizationandControlStructures 15Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).we determine that our desired 13th was or was not a full moon. The variable icon signals the direction ofadjustment. icon=0 1 call flmoon(n,2,jd,frac) Get date off ull moon n. ifrac=nint(24.*(frac+TIMZON)) Convert to hours in correct time zone. if(ifrac.lt.0)then Convert from Julian Days beginning at noon to civil days beginning at midnight. jd=jd-1 ifrac=ifrac+24 endif if(ifrac.gt.12)then jd=jd+1ifrac=ifrac-12 else ifrac=ifrac+12 endifif(jd.eq.jday)then Did we hit our target day? write (*,’(/1x,i2,a,i2,a,i4)’) im,’/’,13,’/’,iyyy write (*,’(1x,a,i2,a)’) ’Full moon ’,ifrac, * ’ hrs after midnight (EST).’ Don’t worry ifyou are unf amiliar with FORTRAN ’s esoteric input/output statements; very few programs in this book do any input/output. goto 2 Part ofthe break-structure, case ofa match. else Didn’t hit it. ic=isign(1,jday-jd) if(ic.eq.-icon) goto 2 Another break, case ofno match. icon=icn=n+ic endif goto 1 2 continue endif enddo 11 enddo 12 END If you are merely curious, there were (or will be) occurrences of a full moon on Friday the 13th (time zone GMT −5) on: 3/13/1903, 10/13/1905, 6/13/1919, 1/13/1922, 11/13/1970,2/13/1987,10/13/2000,9/13/2019, and 8/13/2049. Other “standard” structures. Our advice is to avoid them. Every pro- gramming language has some number of “goodies”that the designer just couldn ’t resist throwing in. They seemed like a good idea at the time. Unfortunately theydon’t stand the testof time! Your program becomes dif ficult to translate into other languages, and dif ficult to read (because rarely used structures are unfamiliar to the reader). You can almost always accomplish the supposed conveniences of these structures in other ways. Try to do so with the above standard structures, which reallyarestandard. If you can ’t, then use straightforward, unstructured, tests and goto’s. This will introduce real (not tame) statement labels, whose very existence will warn the reader to give special thought to the program ’s controlflow. InFORTRAN we consider the ill-advised control structures to be assigned gotoandassignstatements computed gotostatement arithmetic ifstatement 16 Chapter1. PreliminariesSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).About“Advanced Topics” Material set in smaller type, like this, signals an “advanced topic, ”either one outside of the main argument of the chapter, or else one requiring of you more than the usual assumedmathematical background, or else (in a few cases) a discussion that is more speculative or analgorithm that is less well-tested. Nothing important will be lost if you skip the advancedtopics on a first reading of the book. Youmayhavenoticedthat,byitsloopingoverthemonthsandyears,theprogram badluk avoids using any algorithm for converting a Julian Day Number back into a calendar date. Aroutine for doing just this is not very interesting structurally, but it is occasionally useful: SUBROUTINE caldat(julian,mm,id,iyyy) INTEGER id,iyyy,julian,mm,IGREG PARAMETER (IGREG=2299161) I n v e r s eo ft h ef u n c t i o n julday given above. Here julian is input as a Julian Day Number, and the routine outputs mm,id,a n d iyyy as the month, day, and year on which the specified Julian Day started at noon. INTEGER ja,jalpha,jb,jc,jd,jeif(julian.ge.IGREG)then Cross-over to Gregorian Calendar produces this correction. jalpha=int(((julian-1867216)-0.25d0)/36524.25d0) ja=julian+1+jalpha-int(0.25d0*jalpha) else if(julian.lt.0)then Make day number positive by adding in- teger number ofJulian centuries, then subtract them off at the end.ja=julian+36525*(1-julian/36525) else ja=julian endifjb=ja+1524 jc=int(6680.0d0+((jb-2439870)-122.1d0)/365.25d0) jd=365*jc+int(0.25d0*jc)je=int((jb-jd)/30.6001d0) id=jb-jd-int(30.6001d0*je) mm=je-1if(mm.gt.12)mm=mm-12iyyy=jc-4715 if(mm.gt.2)iyyy=iyyy-1 if(iyyy.le.0)iyyy=iyyy-1if(julian.lt.0)iyyy=iyyy-100*(1-julian/36525)return END (Foradditionalcalendricalalgorithms,applicabletovarioushistoricalcalendars,see [8].) Some Habits and Assumed ANSI Extensions Mentioning a few of our programming habits here will make it easier for you to read the programs in this book. We habitually use mandnto refer to the logical dimensions of a matrix, mpandnpto refer to the physical dimensions. (These importantconcepts are detailed in §2.0 and Figure 2.0.1.) Often, when a subroutine or procedure is to be passed some integer n,i t needs an internally preset value for the largest possible value that will be passed. Wehabituallycallthis NMAX,andsetitina PARAMETER statement. Whenwesayinacomment, “largestvalueof n,”wedonotmeantoimply that the program will fail algorithmically for larger values, but only that NMAXmust be altered. A number represented by TINY, usually a parameter, is supposed to be much smaller than any number of interest to you, but not so small that it underflows. Its use is usually prosaic, to prevent divide checks in some circumstances. 1.1ProgramOrganizationandControlStructures 17Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Asamatteroftypography,theprinted FORTRAN programsinthisbook,iftyped into a computerexactly as written, would violate the FORTRAN-77 standard in a few trivialways. Theanomalies,whichare notpresentinthemachine-readableprogram distributions, are as follows: As already discussed, we use enddofollowed by the statement label instead of continue preceded by the label. Standard FORTRAN readsnomorethan72charactersona lineandignores input from column 73 onward. Longer statements are broken up onto “continuation lines. ”In the printed programs in this book, some lines contain more than 72 characters. When the break to a continuation lineis not shown explicitly, it should be inserted when you type the program into a computer. Instandard FORTRAN,columns1through6oneachlineareusedvariously for (i) statement labels, (ii) signaling a comment line, and (iii) signaling a continuation line. We simplify the format slightly: To the left of the “program left margin, ”an integer is a statement label (not a “tame label ” asdescribedabove),anasterisk( *)indicatesacontinuationline,anda “C” indicatesacommentline. Commentlinesshowninthis wayaregenerally either USESstatements(see §1.0),orelse “commented-outprogramlines ” that are separately explained in each instance. A small number of routines in this book require the use of functions that act bitwise on integers, e.g., bitwise “and”or“exclusive or ”. Unfortunately, although thesefunctionsare availableinvirtuallyall modern FORTRAN implementations,they are not a part of the FORTRAN-77 standard. Even more unfortunate is the fact that there are two different naming conventions in widespread use. We use the names iand(i,j) ,ior(i,j) ,not(i),ieor(i,j) , and ishft(i,j) , forand,or,not, exclusive-or , andleft-shift, respectively, as well as the subroutines ibset(i,j) , ibclr(i,j) ,andthelogicalfunction btest(i,j) forbit-set,bit-clear,andbit-test. Some (mainly UNIX) FORTRAN compilers use a different set of names, with the following correspondences: Us... Them ... iand(i,j) =and(i,j) ior(i,j) =or(i,j) not(i) =not(i) ieor(i,j) =xor(i,j) ishft(i,j) =lshft(i,j) ibset(i,j) =bis(j,i) Notereversedarguments! ibclr(i,j) =bic(j,i) Ditto! btest(i,j) =bit(j,i) Ditto! If you are one of “Them,”you can either modify the small number of programs affected (e.g., by inserting FORTRAN statement function de finitions at the beginning of the routines), or else link to an object file into which you have compiled the trivial functions that de fine“our”names in terms of “yours,”as in the above table. Standards really are important! Hexadecimal constants, for which there is no standard notation in FORTRAN compilers,occuratthreeplacesinChapter7: aprogramfragmentattheendof §7.1, 18 Chapter1. PreliminariesSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).and routines psdesandran4in§7.5. We use a notation like Z’3F800000’ , which is consistent with the new FORTRAN-90 standard, but you may need to change this to, e.g., x’3f800000’ ,’3F800000’X ,o re v e n 16#3F800000 . Inextremis,youcan convertthe hexvaluesto decimalintegers;but notethat most compilerswill require anegativedecimalintegeras thevalueofa hexconstantwith its high-orderbitset. Asalreadymentionedin §1.0,thenotation a(1:m),inprogramcommentsandin thetext,denotesthearrayelementrange a(1),a(2),...,a(m). Likewise,notations likeb(2:7)orc(1:m,1:n) aretobeinterpretedasdenotingrangesofarrayindices. CITED REFERENCES AND FURTHER READING: Kernighan, B.W. 1978, The Elements of Programming Style (New York: McGraw-Hill). [1] Yourdon,E.1975, TechniquesofProgramStructureandDesign (EnglewoodCliffs,NJ:Prentice- Hall). [2] Meissner,L.P.andOrganick,E.I.1980, Fortran77FeaturingStructuredProgramming (Reading, MA: Addison-Wesley). [3] Hoare, C.A.R. 1981, Communications of the ACM , vol. 24, pp. 75–83. Wirth, N. 1983, Programming in Modula-2 , 3rd ed. (New York: Springer-Verlag). [4] Stroustrup, B. 1986, The C++Programming Language (Reading, MA: Addison-Wesley). [5] Borland International, Inc. 1989, Turbo Pascal 5.5 Object-Oriented Programming Guide (Scotts Valley, CA: Borland International). [6] Meeus, J. 1982, Astronomical Formulae for Calculators , 2nd ed., revised and enlarged (Rich- mond, VA: Willmann-Bell). [7] Hatcher,D.A.1984, QuarterlyJournaloftheRoyalAstronomicalSociety ,vol.25,pp.53–55;see alsoop. cit.1985, vol. 26, pp. 151–155, and 1986, vol. 27, pp. 506–507. [8] 1.2 Error, Accuracy, and Stability Althoughweassumenopriortrainingofthereaderinformalnumericalanalysis, we will need to presume a common understanding of a few key concepts. We will define these brie fly in this section. Computers store numbers not with in finite precision but rather in some ap- proximation that can be packed into a fixed number of bits(binary digits) or bytes (groups of 8 bits). Almost all computers allow the programmer a choice among several different such representations ordata types . Data types can differ in the number of bits utilized (the wordlength ), but also in the more fundamental respect of whether the stored number is represented in fixed-point (also called integer)o r floating-point (also called real) format. A number in integer representation is exact. Arithmetic between numbers in integerrepresentationisalsoexact,withtheprovisosthat(i)theanswerisnotoutside the range of (usually, signed) integers that can be represented, and (ii) that division is interpretedas producinganintegerresult,throwingaway anyintegerremainder.