/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:58 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sinta.h" #include /* PARAMETER translations */ #define NSIA (NSIB + 1) #define NSIB (NSRA + 1) #define NSRA (NSRB + 1) #define NSRB 1 /* end of PARAMETER translations */ /* COMMON translations */ struct t_sintnc { float ainit, binit, fncval, s, tp, fer, fer1, relobt, tps, xj, xjp; long int fea, fea1, kdim, inc, inc2, istop[2][2], jprint, iprint, kk, kmaxf, ndim, nfindx, nfmax, nfmaxm, reltol, reverm, revers, wherem; LOGICAL32 needh; } sintnc; struct t_sintc { double acum, pacum, result[2]; float aacum, local[4], abscis, ta, delta, delmin, diff, discx[2], end[2], errina, errinb, fat[2], fsave, funct[24], f2, paacum, pf1, pf2, phisum, phtsum, px, space[6], step[2], start[2], sum, t, tasave, tb, tend, worry[2], x1, x2, x, f1, count, xt[17], ft[17], phi[34], absdif, edue2a, edue2b, ep, epnoiz, epsmax, epso, epsr, epss, errat[2], errc, errf, errt[2], esold, extra, pepsmn, releps, rep, rndc, tlen, xjump, erri, err, epsmin, eps, re, reprod; long int discf, dischk, endpts, inew, iold, ip, ixkdim, j, j1, j1old, j2, j2old, kmax, kmin, l, lendt, nfjump, nsubsv, nxkdim, taloc, where2, i, k, kaimt, nsub, part, search, where, nfeval; LOGICAL32 did1, fail, fats[2], fsaved, havdif, iend, init, roundf, xcdobt[2], pad[7]; } sintc; struct t_sintec { float emeps, eepsm8, edelm2, edelm3, esqeps, ersqep, ersqe6, eminf, esmall, enzer, edelm1, eninf; } sintec; /* end of COMMON translations */ void /*FUNCTION*/ sinta( float *answer, float work[], long iflag[]) { LOGICAL32 neweps; long int _d_l, _d_m, _do0, _do1, jumps, konvrg, throw, _i, _r; static long int endphi[2]; float alpha, beta, smin, tp1, zl1; static long aimfor = 6; static float errcf[3]={1.579e0,5.233e-2,8.501e-3}; static long fstore[7]={202,304,508,916,1217,1417,100}; static float fudge[5]={1.0e-8,1.0e-6,1.0e-5,1.0e-4,3.0e-4}; static float gamma[7]={1.42e0,2.28e0,4.62e0,19.00e0,75.00e0,75.00e0, 75.00e0}; static long ih[7]={14,6,16,2,18,8,20}; static long korect[6]={41,81,161,321,55890063,56395705}; static long maxk = 8; static long nsubmx = 4; static float recorr[3]={1.0e-3,0.0e0,1.0e-1}; static float relook = 5.0e-3; static float xeps = 0.02e0; static float xjumps = 0.5e0; static float xstep[7]={0.5e0,0.875e0,1.0e0,2.25e0,5.0e0,15.0e0, 50.0e0}; static float z[4]={1.0e0,1.0e0,1.75e0,6.0e0}; static float p[305]={-.11111111111111111111e00,.22540333075851662296e00, .55555555555555555556e00,.647209421402969791e-02,-.928968790944433705e-02, .39508731291979716579e-01,.10465622602646726519e00,.56575625065319744200e00, .40139741477596222291e00,.5223046896961622e-04,.17121030961750000e-03, -.724830016153892898e-03,-.7017801099209042e-04,.61680367872449777899e-02, .17001719629940260339e-01,.11154076712774300110e00,.92927195315124537686e-01, .37889705326277359705e00,.17151190913639138079e00,.77661331357103311837e00, .21915685840158749640e00,.682166534792e-08,.12667409859336e-06, .59565976367837165e-05,.1392330106826e-07,-.6629407564902392e-04, -.704395804282302e-06,-.34518205339241e-07,-.814486910996e-08, .90187503233240234038e-03,.25447807915618744154e-02,.18468850446259893130e-01, .16446049854387810934e-01,.70345142570259943330e-01,.35957103307129322097e-01, .16327406183113126449e00,.56979509494123357412e-01,.29750379350847292139e00, .76879620499003531043e-01,.46868025635562437602e00,.93627109981264473617e-01, .66886460674202316691e00,.10566989358023480974e00,.88751105686681337425e00, .11195687302095345688e00,.371583e-15,.21237877e-12,.10522629388435e-08, .1748029e-14,.3475718983017160e-06,.90312761725e-11,.12558916e-13, .54591e-15,-.72338395508691963e-05,-.169699579757977e-07,-.854363907155e-10, -.12281300930e-11,-.462334825e-13,-.42244055e-14,-.88501e-15, -.40904e-15,.12711187964238806027e-03,.36322148184553065969e-03, .27937406277780409196e-02,.25790497946856882724e-02,.11315242452570520059e-01, .61155068221172463397e-02,.27817125251418203419e-01,.10498246909621321898e-01, .53657141626597094849e-01,.15406750466559497802e-01,.89628843042995707499e-01, .20594233915912711149e-01,.13609206180630952284e00,.25869679327214746911e-01, .19305946804978238813e00,.31073551111687964880e-01,.26024395564730524132e00, .36064432780782572640e-01,.33709033997521940454e00,.40715510116944318934e-01, .42280428994795418516e00,.44914531653632197414e-01,.51638197305415897244e00, .48564330406673198716e-01,.61664067580126965307e00,.51583253952048458777e-01, .72225017797817568492e00,.53905499335266063927e-01,.83176474844779253501e00, .55481404356559363988e-01,.94365568695340721002e00,.56277699831254301273e-01, .1041098e-15,.249472054598e-10,.55e-20,.290412475995385e-07,.367282126e-13, .5568e-18,-.871176477376972025e-06,-.8147324267441e-09,-.8830920337e-12, -.18018239e-14,-.70528e-17,-.506e-19,.17569645108401419961e-04, .50536095207862517625e-04,.40120032808931675009e-03,.37774664632698466027e-03, .16833646815926074696e-02,.93836984854238150079e-03,.42758953015928114900e-02, .16811428654214699063e-02,.85042788218938676006e-02,.25687649437940203731e-02, .14628500401479628890e-01,.35728927835172996494e-02,.22858485360294285840e-01, .46710503721143217474e-02,.33362148441583432910e-01,.58434498758356395076e-02, .46269993574238863589e-01,.70724899954335554680e-02,.61679602220407116350e-01, .83428387539681577056e-02,.79659974529987579270e-01,.96411777297025366953e-02, .10025510022305996335e00,.10955733387837901648e-01,.12348658551529473026e00, .12275830560082770087e-01,.14935550523164972024e00,.13591571009765546790e-01, .17784374563501959262e00,.14893641664815182035e-01,.20891506620015163857e00, .16173218729577719942e-01,.24251603361948636206e00,.17421930159464173747e-01, .27857691462990108452e00,.18631848256138790186e-01,.31701256890892077191e00, .19795495048097499488e-01,.35772335749024048622e00,.20905851445812023852e-01, .40059606975775710702e00,.21956366305317824939e-01,.44550486736806745112e00, .22940964229387748761e-01,.49231224246628339785e00,.23854052106038540080e-01, .54086998801016766712e00,.24690524744487676909e-01,.59102017877011132759e00, .25445769965464765813e-01,.64259616216846784762e00,.26115673376706097680e-01, .69542355844328595666e00,.26696622927450359906e-01,.74932126969651682339e00, .27185513229624791819e-01,.80410249728889984607e00,.27579749566481873035e-01, .85957576684743982540e00,.27877251476613701609e-01,.91554595991628911629e00, .28076455793817246607e-01,.97181535105025430566e00,.28176319033016602131e-01, .3326e-18,.114094770478e-11,.2952436056970351e-08,.51608328e-15, -.110177219650597323e-06,-.58656987416475e-10,-.23340340645e-13, -.1248950e-16,.24036202515353807630e-05,.69379364324108267170e-05, .56003792945624240417e-04,.53275293669780613125e-04,.23950907556795267013e-03, .13575491094922871973e-03,.61966197497641806982e-03,.24921240048299729402e-03, .12543855319048853002e-02,.38974528447328229322e-03,.21946455040427254399e-02, .55429531493037471492e-03,.34858540851097261500e-02,.74028280424450333046e-03, .51684971993789994803e-02,.94536151685852538246e-03,.72786557172113846706e-02, .11674841174299594077e-02,.98486295992298408193e-02,.14049079956551446427e-02, .12907472045965932809e-01,.16561127281544526052e-02,.16481342421367271240e-01, .19197129710138724125e-02,.20593718329137316189e-01,.21944069253638388388e-02, .25265540247597332240e-01,.24789582266575679307e-02,.30515340497540768229e-01, .27721957645934509940e-02,.36359378430187867480e-01,.30730184347025783234e-02, .42811783890139037259e-01,.33803979910869203823e-02,.49884702478705123440e-01, .36933779170256508183e-02,.57588434808916940190e-01,.40110687240750233989e-02, .65931563842274211999e-01,.43326409680929828545e-02,.74921067092924347640e-01, .46573172997568547773e-02,.84562412844234959360e-01,.49843645647655386012e-02, .94859641186738404810e-01,.53130866051870565663e-02,.10581543166444097714e00, .56428181013844441585e-02,.11743115975265809315e00,.59729195655081658049e-02, .12970694445188609414e00,.63027734490857587172e-02,.14264168911376784347e00, .66317812429018878941e-02,.15623311732729139895e00,.69593614093904229394e-02, .17047780536259859981e00,.72849479805538070639e-02,.18537121234486258656e00, .76079896657190565832e-02,.20090770903915859819e00,.79279493342948491103e-02, .21708060588171698360e00,.82443037630328680306e-02,.23388218069623990928e00, .85565435613076896192e-02,.25130370638306339718e00,.88641732094824942641e-02, .26933547875781873867e00,.91667111635607884067e-02,.28796684463774796540e00, .94636899938300652943e-02,.30718623022088529711e00,.97546565363174114611e-02, .32698116976958152079e00,.10039172044056840798e-01,.34733833458998250389e00, .10316812330947621682e-01,.36824356228880576959e00,.10587167904885197931e-01, .38968188628481359983e00,.10849844089337314099e-01,.41163756555233745857e00, .11104461134006926537e-01,.43409411457634557737e00,.11350654315980596602e-01, .45703433350168850951e00,.11588074033043952568e-01,.48044033846254297801e00, .11816385890830235763e-01,.50429359208123853983e00,.12035270785279562630e-01, .52857493412834112307e00,.12244424981611985899e-01,.55326461233797152625e00, .12443560190714035263e-01,.57834231337383669993e00,.12632403643542078765e-01, .60378719394238406082e00,.12810698163877361967e-01,.62957791204992176986e00, .12978202239537399286e-01,.65569265840056197721e00,.13134690091960152836e-01, .68210918793152331682e00,.13279951743930530650e-01,.70880485148175331803e00, .13413793085110098513e-01,.73575662758907323806e00,.13536035934956213614e-01, .76294115441017027278e00,.13646518102571291428e-01,.79033476175681880523e00, .13745093443001896632e-01,.81791350324074780175e00,.13831631909506428676e-01, .84565318851862189130e00,.13906019601325461264e-01,.87352941562769803314e00, .13968158806516938516e-01,.90151760340188079791e00,.14017968039456608810e-01, .92959302395714482093e00,.14055382072649964277e-01,.95773083523463639678e00, .14080351962553661325e-01,.98590611358921753738e00,.14092845069160408355e-01}; static int _aini = 1; float *const alocal = (float*)sintc.local; float *const blocal = (float*)((float*)sintc.local + 1); float *const fata = (float*)sintc.fat; LOGICAL32 *const fatas = (LOGICAL32*)sintc.fats; float *const fatb = (float*)((float*)sintc.fat + 1); LOGICAL32 *const fatbs = (LOGICAL32*)((LOGICAL32*)sintc.fats + 1); long int *const ic1 = (long*)sintnc.istop; long int *const ic2 = (long*)((long*)sintnc.istop + 2); float (*const phis)[17] = (float(*)[17])sintc.phi; float *const phit = (float*)((float*)sintc.phi + 17); /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Discx = &sintc.discx[0] - 1; float *const End = &sintc.end[0] - 1; long *const Endphi = &endphi[0] - 1; float *const Errat = &sintc.errat[0] - 1; float *const Errcf = &errcf[0] - 1; float *const Errt = &sintc.errt[0] - 1; float *const Fat = &sintc.fat[0] - 1; LOGICAL32 *const Fats = &sintc.fats[0] - 1; long *const Fstore = &fstore[0] - 1; float *const Ft = &sintc.ft[0] - 1; float *const Fudge = &fudge[0] - 1; float *const Funct = &sintc.funct[0] - 1; float *const Gamma = &gamma[0] - 1; long *const Iflag = &iflag[0] - 1; long *const Ih = &ih[0] - 1; long *const Korect = &korect[0] - 1; float *const Local = &sintc.local[0] - 1; float *const P = &p[0] - 1; float *const Phi = &sintc.phi[0] - 1; float *const Phit = &phit[0] - 1; float *const Recorr = &recorr[0] - 1; double *const Result = &sintc.result[0] - 1; float *const Start = &sintc.start[0] - 1; float *const Step = &sintc.step[0] - 1; float *const Work = &work[0] - 1; float *const Worry = &sintc.worry[0] - 1; LOGICAL32 *const Xcdobt = &sintc.xcdobt[0] - 1; float *const Xstep = &xstep[0] - 1; float *const Xt = &sintc.xt[0] - 1; float *const Z = &z[0] - 1; /* end of OFFSET VECTORS */ if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ Endphi[1] = 1; _aini = 0; } /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 2009-11/03 SINTA Krogh Set initial value for a few variables. *>> 2009-10-17 SINTA Krogh AACUM now set with ACUM on coalescing nodes. *>> 2008-02-23 SINTA Krogh/Snyder Misc. Changes *>> 2001-09-11 SINTA Krogh Fixed so really small intervals not picked. *>> 1996-04-27 SINTA Krogh Changes to use .C. and C%%. *>> 1996-03-31 SINTA Krogh Removed unused variable in common. *>> 1996-03-30 SINTA Krogh Change specific intrinsics to generics. *>> 1995-11-28 SINTA Krogh Converted from SFTRAN to Fortran 77. *>> 1995-10-24 SINTA Krogh Removed blanks in numbers for C conversion. *>> 1994-11-14 SINTA Krogh Declared all vars. *>> 1994-10-19 SINTA Krogh Changes to use M77CON *>> 1994-09-01 SINTA Snyder don't give round-off msg if err << eps *>> 1994-08-22 SINTA Krogh -- Modified data for C conversion. *>> 1994-08-19 SINTA Snyder fix divide by 0 if ERR=0 and don't believe. *>> 1994-08-15 SINTA Snyder convert MAX(A,B,C) to MAX(A,MAX(B,C)) *>> 1994-07-07 SINTA Snyder set up for CHGTYP. *>> 1994-06-29 SINTA Snyder suppress discontinuity msg if one jump. *>> 1994-06-22 SINTA Snyder don't use RE when reducing stepsize. *>> 1993-05-26 SINTA Krogh -- Added data for C conversion. Needs edit. *>> 1993-05-18 SINTA Krogh -- Changed "END" to "END PROGRAM" *>> 1992-03-03 SINTA Added call to SINTO for an error message. *>> 1991-09-20 SINTA Krogh converted '(1)' dimensioning to '(*)'. *>> 1989-11-08 SINTA Snyder Revise coalesced abscissa test for CRAY *>> 1988-05-17 SINTA Snyder Initial code. *--S replaces "?": ?INT, ?INTA, ?INTC, ?INTDL, ?INTDU, ?INTEC *--& ?INTF, ?INTM, ?INTMA, ?INTNC, ?INTNS, ?INTO, ?INTSM, ?INT1 * * THIS IS THE INTEGRATION PROGRAM FOR SINT. * THE RESULT IS OBTAINED USING QUADRATURE FORMULAE DUE TO * T. N. L. PATTERSON, MATHEMATICS OF COMPUTATION, VOLUME 22, * PAGES 847-856, 1968. * * ***** WARNING *********************************************** * * THE RELIABILITY AND EFFICIENCY OF THIS PROGRAM ARE STRONGLY * INFLUENCED BY DISCONTINUITIES IN THE FUNCTION OR IT'S * DERIVATIVES, INTEGRABLE SINGULARITIES ON THE REAL AXIS, * NON-INTEGRABLE SINGULARITIES NEAR THE REGION, AND POLES * NOT ON THE REAL AXIS. THE EFFICIENCY AND RELIABILITY * OF INTEGRATING SUCH FUNCTIONS MAY BE GREATLY IMPROVED BY MANUALLY * SUBDIVIDING THE INTERVAL OF INTEGRATION AT THE DISCONTINUITY, * SINGULARITY OR REAL PART OF THE POLE, AND SUMMING THE ANSWERS. * A CHANGE OF VARIABLE TO ELIMINATE OR REDUCE THE STRENGTH OF THE * SINGULARITY WILL SIGNIFICANTLY IMPROVE PERFORMANCE. * * ***** FORMAL ARGUMENTS ************************************* * * ANSWER IS THE COMPUTED INTEGRAL IF IFLAG .LT. 0. IF IFLAG = 6, * ANSWER IS THE APPROXIMATE LOCATION OF A SINGULARITY WHICH * APPEARS NOT TO BE INTEGRABLE. IF IFLAG = 0, ANSWER IS * USED TO PASS A FUNCTION VALUE FROM THE USER PROGRAM TO SINTA. */ /* WORK CONTAINS THE ESTIMATED ABSOLUTE ERROR IF IFLAG .LT. 0. * IF IFLAG = 0, WORK IS USED TO PASS AN ABSCISSA TO THE * USER PROGRAM (FOR REVERSE COMMUNICATION). */ /* IFLAG IS USED TO CONTROL REVERSE COMMUNICATION AND TO RETURN * STATUS INFORMATION TO THE USER. VALUES ARE * 0 - INTEGRAND VALUE NEEDED. COMPUTE ANSWER=INTEGRAND * EVALUATED AT WORK(1), ..., WORK(NDIM). * -1 - NORMAL TERMINATION, EITHER THE ABSOLUTE OR THE RELATIVE * ERROR CRITERIA ARE SATISFIED. * -2 - NORMAL TERMINATION, NEITHER THE ABSOLUTE NOR RELATIVE * ERROR CRITERIA ARE SATISFIED, BUT THE RELATIVE TO THE * OBTAINABLE ERROR CRITERION IS SATISFIED. * -3 - NORMAL TERMINATION, BUT NONE OF THE ERROR CRITERIA ARE * SATISFIED. * 4 - BAD IOPT VALUE (SEE SINT1 OR SINTM). * 5 - TOO MANY FUNCTION VALUES NEEDED. * KDIM+5 - APPARENT NON-INTEGRABLE SINGULARITY NEAR ANSWER. * KDIM IS CURRENT DIMENSION - 1 IN 1 DIMENSIONAL CASE. */ /* ***** EXTERNAL REFERENCES *********************************** * * SINTDL TO FORM AND ANALYZE DIFFERENCE LINES. * * SINTDU TO UPDATE DIFFERENCE LINES. * * SINTF PROVIDE THE INTEGRAND. SINTF IS EXPLAINED IN * EXPLAINED IN SINT1 AND SINTM. * * SINTNS TO INCREASE AND REDUCE NSUB, THE DEGREE OF THE INDEPENDENT * VARIABLE TRANSFORMATION. * * SINTO USED FOR PRINTING OUTPUT. * * SINTSM USED TO CALCULATE THE MINIMUM STEPSIZE. */ /* ***** PARAMETERS ***************************************** * * NSIA AS FOR NSRA, BUT NSUB IS TO BE INCREASED. * NSIB AS FOR NSRB, BUT NSUB IS TO BE INCREASED. * NSRA IS AN ARGUMENT OF SINTNS THAT MEANS REDUCE NSUB, BUT ONLY * ALOCAL NEEDS TO BE RECALCULATED. * NSRB AS FOR NSRA, BUT BOTH ALOCAL AND BLOCAL MUST BE RECALCULATED. * */ /* ***** INTERNAL AND COMMON VARIABLES ********************** * * AACUM IS THE INTEGRAL OF THE ABSOLUTE VALUE OF THE INTEGRAND OVER * CURRENT PANEL. * ABSCIS IS THE ABSCISSA TO BE TRANSFORMED TO USER COORDINATES. * SEE NSUB. * ABSDIF IS THE ABSOLUTE VALUE OF ONE HALF OF THE LENGTH OF THE PANEL * (SMALLEST SUBDIVISION) CURRENTLY BEING INTEGRATED. * ACUM IS THE ESTIMATE OF THE INTEGRAL OVER THE CURRENT PANEL. * AIMFOR IS AN INDEX OF THE FORMULA TO AIM FOR AFTER A PANEL HAS BEEN * SUCCESSFULLY INTEGRATED. THE CURRENT VALUE OF AIMFOR * SPECIFIES A 63-POINT FORMULA. */ /* AINIT IS THE INITIAL LOWER REGION BOUNDARY. * ALOCAL IS THE 'LEFT' BOUNDARY OF THE PANEL BEING INTEGRATED. */ /* ALPHA IS USED DURING DETERMINATION OF THE NEXT STEP SIZE AS AN * INDICATOR OF THE PERFORMANCE OF THE INTEGRATION OF THE * PREVIOUS STEP, WITH RESPECT TO THE COMMITTED AND PERMITTED * ERROR. */ /* BETA IS USED DURING COMPUTATION OF ABSCISSAE USED BY THE * QUADRATURE FORMULA AND DURING COMPUTATION OF THE NEXT STEP. */ /* BINIT IS THE INITIAL UPPER REGION BOUNDARY. * BLOCAL IS THE 'RIGHT' BOUNDARY OF THE PANEL CURRENTLY BEING * INTEGRATED. */ /* COUNT IS USED DURING A SEARCH TO DISTINGUISH INTEGRABLE FROM * NON-INTEGRABLE SINGULARITIES. * DELMIN IS THE MINIMUM STEPSIZE ALLOWED, LEST THE ABSCISSAE COALESCE. * DELTA IS THE ABSOLUTE VALUE OF THE CURRENT STEP SIZE. * DID1 IS INITIALLY FALSE, AND SET TRUE WHEN ONE PANEL HAS * BEEN SUCCESSFULLY INTEGRATED. * DIFF IS ONE HALF OF THE SIGNED CURRENT STEP SIZE. * DISCF IS NON-ZERO WHEN A DISCONTINUITY IS DETECTED. IF DISCF IS * POSITIVE AND IPRINT IS NOT ZERO AFTER SUCCESSFULLY * INTEGRATING A PANEL IN PART = 1, THE BOUNDARIES FOR THE * PANEL ARE PRINTED. IN ADDITION, IF ABS(DISCF) IS 2 AFTER A * PANEL IS INTEGRATED, THE REST OF THE PART IS USED FOR THE * NEXT PANEL. * DISCHK IS SET TO MAX(DISCHK,1) WHEN A JUMP IS SEEN IN THE DIFFERENCE * LINE, BUT THE NEXT FORMULA IS TRIED, OR WHEN A DISCONTINUITY * IS DETECTED. IF THERE ARE JUMPS IN THE DIFFERENCE LINE WHICH * SUBSEQUENTLY VANISH, DISCHK IS SET TO 2. WHEN A PANEL IS * ACCEPTED, DISCHK IS DECREMENTED AS FAR AS ZERO. DISCHK IS * SET TO -1 WHEN NSUB IS CHANGED. DISCHK IS SET TO -2 WHEN A * CAUTION POINT IS ENCOUNTERED. IF DISCHK IS POSITIVE, THE * DIFFERENCE LINES MUST BE FORMED EVERY QUADRATURE STEP, AND * THE END POINTS MUST BE ADDED BEFORE THE ANSWER CAN BE * ACCEPTED. IF DISCHK IS -1, THE DIFFERENCE LINE MUST BE * FORMED BEFORE THE ANSWER CAN BE ACCEPTED. IF DISCHK IS -2 * THE DIFFERENCE LINES MUST BE FORMED AND THE END POINTS ADDED * BEFORE THE ANSWER CAN BE ACCEPTED. * DISCX IS THE APPROXIMATE LOCATION OF A DISCONTINUITY. * EDUE2A, EDUE2B ARE THE ERRORS IN THE INTEGRAL DUE TO IMPRECISE LIMITS * A AND B RESPECTIVELY. * END STORES THE 'RIGHT' BOUNDARIES OF MAJOR SUBDIVISIONS. * ENDPHI CONTAINS THE END OF THE DIFFERENCE LINE ON WHICH THERE CAN BE * A JUMP. ENDPHI(1) = 1 BECAUSE THERE CAN ONLY BE A JUMP IN PHI * AT 1. ENDPHI(2) = LENDT BECAUSE THERE CAN ONLY BE A JUMP IN * PHIT AT LENDT. ENDPHI IS USED ONLY WHEN CHECKING WHETHER THE * ERROR IN THE FUNCTION AT THE END OF A PANEL IS SO LARGE THAT * A JUMP IS NOT NECESSARILY BELIEVABLE. */ /* ENDPTS IS USED TO CONTROL ADDITION OF THE END POINTS TO THE * DIFFERENCE LINES: * 1. FUNCTION NEED NOT BE EVALUATED AT A. * 2. FUNCTION NEED NOT BE EVALUATED AT B. * 3. FUNCTION NEED NOT BE EVALUATED AT EITHER A AND B. * EP IS THE PREVIOUS VALUE OF ERRC. * EPNOIZ IS THE MAXIMUM ERROR THRESHOLD FOR NOISE DETECTION. * EPS IS THE ERROR ALLOWED ON THE CURRENT PANEL. * EPSMAX IS THE MOST STRINGENT TOLERANCE ALLOWED ON THE CURRENT PANEL. * SEE XEPS. * EPSMIN IS AN ESTIMATE OF THE SMALLEST ERROR POSSIBLE ON THE CURRENT * PANEL. * EPSO IS THE ORIGINAL ERROR TOLERANCE. * EPSR IS THE REMAINING ERROR TOLERANCE. * EPSS IS A STORAGE CELL FOR THE MAXIMUM OF EPS AND EPSMAX. EPS IS * CHANGED AFTER EPSS IS COMPUTED. * ERR IS THE EXTRAPOLATED ERROR COMMITTED DURING THE CURRENT * INTEGRATION STEP. * ERRAT ERRAT(I) CONTAINS THE ERROR IN FAT(I) IF FATAS(I), I=1 OR 2. * ERRC IS THE ERROR COMMITTED DURING THE CURRENT INTEGRATION STEP, * FOR PURPOSES OF DETERMINING THE NEXT STEPSIZE. * ERRCF IS USED TO COMPUTE A MINIMUM VALUE FOR THE ERROR * WHEN K.LE.4, USING THE DIFFERENCE LINES. */ /* ERRF IS THE ESTIMATED ERROR IN THE INTEGRAND. WHEN KDIM=1, ERRF * IS OBTAINED FROM WORK(FEA) (IF FEA .NE. 0). WHEN KDIM.GT.1 * ERRF IS THE ESTIMATED ERROR IN THE INNER INTEGRAL. * ERRI IS THE ERROR ESTIMATED BY THE DIFFERENCE BETWEEN TWO * SUCCESSIVE INTEGRATION FORMULAE APPLIED TO THE CURRENT PANEL. * ERRINA, ERRINB ARE THE ERRORS THE USER DECLARES BY WAY OF AN OPTION * TO BE PRESENT IN THE LIMITS A AND B RESPECTIVELY. * ERRT IS THE TOTAL ERROR COMMITTED ON A MAJOR SUBDIVISION. * ESOLD IS THE PREVIOUS ERROR ESTIMATED ACROSS A JUMP BY THE SEARCH. * EXTRA IS ANY EXTRA AMOUNT ADDED TO THE ERROR BECAUSE OF OUR * SUSPICIOUS NATURE. EXTRA IS BASED UPON HIGH-ORDER * DIFFERENCES IF THERE ARE NO JUMPS. * FAIL INDICATES TO THE STEPSIZE SELECTION ROUTINE THAT INTEGRATION * WAS NOT SUCCESSFUL. FAIL IS ALSO SET IF A JUMP IN THE * DIFFERENCE LINE VANISHES WHEN ACCURATE LOCATION IS ATTEMPTED. * FAT IS A VECTOR EQUIVALENT TO FATA,FATB. * FATA SAVES THE FUNCTION VALUE AT ALOCAL. */ /* FATAS IS TRUE WHEN FATA AND ERRAT(1) CONTAIN VALID DATA. */ /* FATB SAVES THE FUNCTION VALUE AT BLOCAL. */ /* FATBS IS TRUE WHEN FATB AND ERRAT(2) CONTAIN VALID DATA. */ /* FATS IS A VECTOR EQUIVALENT TO FATAS,FATBS. * FEA IS THE INDEX IN WORK OF THE ABSOLUTE ERROR COMMITTED * COMPUTING F. WORK(FEA) MAY BE CHANGED DURING THE INTEGRATION. * FEA APPLIES ONLY TO THE INNERMOST INTEGRATION OF A MULTIPLE * QUADRATURE. * FER IS THE RELATIVE ERROR COMMITTED COMPUTING F. * FNCVAL IS THE COMPUTED FUNCTION VALUE. * FSAVE IS USED TO STORE THE FUNCTION VALUE AT THE ABSCISSA WHERE * THE INTERVAL IS SUBDIVIDED. * FSAVED IS SET IF FSAVE CONTAINS USEFUL DATA. * FSTORE IS USED TO DETERMINE WHICH FUNCTION VALUES ARE TO BE STORED, * AND WHERE THEY ARE TO BE STORED IN THE FUNCTION TABLE. * FSTORE IS READ-ONLY. THE CONTENT OF EACH ELEMENT OF FSTORE * IS A TWO DIGIT BASE-100 NUMBER, WHERE THE FIRST DIGIT * INDICATES WHERE THE FIRST FUNCTION VALUE COMPUTED DURING THE * CURRENT INTEGRATION STEP IS TO BE STORED, AND THE SECOND DIGIT * INDICATES THE LAST CELL OF THE FUNCTION TABLE TO BE USED FOR * THE CURRENT INTEGRATION STEP. */ /* FT USED TO HOLD, TEMPORARILY, THE FUNCTION VALUES USED. * THESE FUNCTION VALUES ARE USED TO COMPUTE THE DIFFERENCE * LINE. SEE XT. * FUDGE IS USED DURING THE CONVERGENCE TEST TO RELATE THE ESTIMATED * ERROR TO THE ROUND-OFF LEVEL. FUDGE(K-KAIMT) GIVES THE * RELATIVE PRECISION THAT MUST BE OBTAINED IF THE ERROR * ESTIMATE IS NOT TO BE INCREASED. */ /* FUNCT IS USED TO STORE FUNCTION VALUES WHICH MAY BE RE-USED BY * THE PATTERSON FORMULAE, AND ADDITIONAL FUNCTION VALUES TO * COMPUTE THE DIFFERENCE LINE. * FUNCT(1-17) IS STORAGE NEEDED FOR THE PATTERSON METHOD * FUNCT(18-24) IS STORAGE TO RECOVER F(X1) THROUGH F(X15). * F1, F2 CONTAIN NEWLY COMPUTED FUNCTION VALUES DURING THE QUADRATURE * STEP. THEY ARE ALSO USED ELSEWHERE FOR TEMPORARY STORAGE. * GAMMA IS USED DURING COMPUTATION OF THE NEXT STEPSIZE. */ /* HAVDIF IS SET TO INDICATE THAT THE DIFFERENCE LINE HAS BEEN COMPUTED * USING 15 FUNCTION VALUES. * I IS USED FREELY AS AN INDEX. * IC1 IC2 ARE THE INDICES OF CONVERGENCE IN THE FORWARD AND BACKWARD * DIFFERENCE LINES. IC1 AND IC2 ARE EQUIVALENT TO ISTOP(1,1) * AND ISTOP(1,2). (SEE BELOW). * IEND IS SET WHEN THE END OF A MAJOR SUBDIVISION IS IN THE CURRENT * PANEL. * IH ARE INDICES USED TO LOCATE ABSCISSAE IN THE TABLE OF * QUADRATURE RULES. IH IS USED DURING COMPUTATION OF THE * DIFFERENCE LINE. */ /* INC INC2 ARE USED DURING ANALYSIS OF THE DIFFERENCE LINE AS INDEX * INCREMENTORS. THEY ARE ALSO USED FREELY AS INDICES. * INEW IS USED DURING THE QUADRATURE STEP TO CONTROL A LOOP. INEW * MAY NOT BE CHANGED BETWEEN QUADRATURE STEPS. * INIT IS SET TO INDICATE THAT THE ORIGINAL INTERVAL IS BEING * INTEGRATED. * INSTOP IS USED TO CONTROL A LOOP DURING ANALYSIS OF THE DIFFERENCE * LINE. * IOLD IS USED DURING THE QUADRATURE STEP TO CONTROL A LOOP. IOLD * MAY NOT BE CHANGED BETWEEN QUADRATURE STEPS. * IP IS USED DURING THE QUADRATURE STEP TO INDEX THE TABLE OF * ABSCISSAE AND WEIGHTS. IP MAY NOT BE CHANGED BETWEEN * QUADRATURE STEPS. * IPRINT CONTROLS THE LEVEL OF DIAGNOSTIC PRINTING. * 0 - NO PRINTING * 1 - MINIMUM PRINTING - ERROR MESSAGES (DEFAULT) * 2 - PANEL BOUNDARIES AND ANSWERS * 3 - ERROR ESTIMATES FOR EACH QUADRATURE FORMULA * 4 - DETAILED OUTPUT (DIFFERENCE LINES, ETC). * ISTOP HOLDS RESULTS OF ANALYSIS OF DIFFERENCE LINE. * ENTRIES IN ISTOP MEAN * ISTOP(1,1) - INDEX OF XT TO WHICH CONVERGENCE EXTENDS FROM * THE ALOCAL END OF THE INTERVAL, * ISTOP(1,2) - INDEX OF XT TO WHICH CONVERGENCE EXTENDS FROM * THE BLOCAL END OF THE INTERVAL * ISTOP(2,1) AND ISTOP(2,2) MAY INDICATE INDICES OF XT WHERE * THE INTERVAL SHOULD BE SUBDIVIDED. * IXKDIM INDEXES THE CELL IN IFLAG IN WHICH THE CURRENT DIMENSION OF A * MULTIDIMENSIONAL INTEGRATION IS STORED. THE DEFAULT IS 1, * CORRESPONDING WITH THE USUAL CELL USED FOR COMMUNICATION WITH * SINTF. THE USER MAY NEED TO CHANGE THIS BY USING OPTION 12 * IF NON-STANDARD DIMENSION CHANGES ARE NECESSARY, AND A CHANGE * OF VARIABLE IS SIMULTANEOUSLY NECESSARY. SEE SINTM AND * SINTMA. * J IS FREELY USED AS AN INDEX. * JUMPS STORES THE QUALITATIVE RESULT OF THE DIFFERENCE LINE ANALYSIS. * 1 - NO JUMPS * 2 - ONE JUMP, ON THE END OF THE LINE * 3 - ONE JUMP, NOT ON THE END * 4 - TWO JUMPS. * J1 J2 ARE USED DURING THE QUADRATURE STEP TO CONTROL STORING INTO * THE FUNCT ARRAY. THEY ARE RECOMPUTED AT EACH STEP. THEY * ARE ALSO USED DURING THE ANALYSIS OF THE DIFFERENCE LINE TO * STORE 'JUMP' INDICES. * J1OLD ARE USED TO REMEMBER J1 AND J2 DURING A SEARCH * J2OLD FOR A JUMP. * K IS THE INDEX FOR THE QUADRATURE FORMULA LAST USED. THE * NUMBER OF POINTS USED IN THE KTH FORMULA IS 2**K - 1. * KAIMT IS THE INDEX OF THE QUADRATURE FORMULA EXPECTED TO CONVERGE * ON THE NEXT STEP. * KDIM IS THE CURRENT DIMENSION - ALWAYS 1 IN 1 DIMENSIONAL CASE. * KK IS A TEMPORARY INTEGER STORAGE. * KMAX IS THE MAXIMUM VALUE ALLOWED FOR K ON THE CURRENT INTEGRATION * STEP. * KMAXF IS THE VALUE TO USE FOR KMAX INITIALLY AND AFTER A FAILURE. * KMIN IS THE MINIMUM QUADRATURE FORMULA INDEX TO ACCEPT IF THERE * ARE NO JUMPS. KMIN IS USUALLY ZERO, BUT IS SET TO K IF THE * STEP SIZE INCREASES. * KONVRG KONVRG IS THE TOTAL NUMBER OF NODES OF CONVERGENCE OBSERVED * IN THE DIFFERENCE LINES. KONVRG=IC1+LENDT-IC2. * KORECT IS USED TO DETERMINE WHICH OF THE STORED FUNCTION VALUES ARE * TO BE USED TO CORRECT ONE-HALF OF THE PREVIOUSLY COMPUTED * INTEGRAL ON THE CURRENT INTERVAL. THE INDICES ARE STORED IN * KORECT AS BASE-20 DIGITS. KORECT IS READ-ONLY, INDEXED BY K. */ /* L IS USED DURING ANALYSIS OF THE DIFFERENCE LINE. * LENDT IS THE CURRENT LENGTH OF THE DIFFERENCE LINE. * LOCAL IS A VECTOR EQUIVALENT TO ALOCAL,BLOCAL. * LOCAL(3),LOCAL(4) ARE ALOCAL,BLOCAL IN USER COORDINATES. * MAXK IS THE MAXIMUM VALUE ALLOWED FOR K. * NEEDH IS SET IF THE HEADER FOR DIAGNOSTIC INFORMATION IS NEEDED. * NEWEPS MEANS EPS MUST BE RECOMPUTED TO ACCOUNT FOR ERROR IN LIMITS */ /* NFEVAL IS THE CURRENT NUMBER OF FUNCTION EVALUATIONS. * NFJUMP IS THE VALUE OF NFEVAL WHEN A JUMP IS FIRST NOTICED. * NFMAX IS THE INDEX OF THE ENTRY IN THE OPTION VECTOR (IOPT) TO * USE TO CONTROL THE MAXIMUM NUMBER OF FUNCTION EVALUATIONS. * NOMOUT IS THE NOMINAL OUTPUT UNIT. * NSUB IS THE POWER OF THE TRANSFORMATION OF THE FORM * X=TA+(T-TA)*((T-TA)/(TB-TA))**(NSUB-1) * WHERE X IS THE USER'S INDEPENDENT VARIABLE AND T IS THE * INDEPENDENT VARIABLE TO USE FOR INTEGRATION. NSUB MAY * ONLY BE A POWER OF 2. NSUB=0 REALLY MEANS NSUB=1, BUT * IS SIMPLER TO USE. * NSUBMX IS THE MAXIMUM VALUE FOR NSUB. * NSUBSV IS THE VALUE TO USE FOR NSUB WHEN RESTARTING * AFTER A SUBDIVISION. * NXKDIM IS THE VALUE OF KDIM BEFORE AN INNER INTEGRAL WAS REQUESTED. * THE INNER INTEGRAL IS NOT NECESSARILY OF DIMENSION KDIM-1. * SEE IXKDIM ABOVE, SINTM AND SINTMA. * P IS THE ARRAY OF ABSCISSAE AND WEIGHTS. DATA FOR P IS INCLUDED * AFTER THE DESCRIPTION OF INTERNAL VARIABLES. */ /* PAACUM IS THE PREVIOUS VALUE OF AACUM. * PACUM IS THE PREVIOUS VALUE OF ACUM. ABS(ACUM-PACUM) IS USED AS AN * INITIAL ERROR ESTIMATE. * PART IS THE INDEX OF THE CURRENT MAJOR SUBDIVISION. START, END, * STEP, ERRT AND RESULT ARE SUBSCRIPTED BY PART. PART BEGINS * SET TO 1, IS CHANGED TO 2 WHEN THE INTERVAL IS SUBDIVIDED, * AND SET BACK TO 1 WHEN PART 2 IS INTEGRATED. WHEN PART 1 * IS FINISHED, THE ENTIRE PROBLEM IS FINISHED. * PEPSMN IS THE PREVIOUS VALUE OF EPSMIN. PEPSMN MAY NOT BE CHANGED * BETWEEN QUADRATURE STEPS. * PERR IS THE ERROR COMMITTED ON THE LAST SUCCESSFULLY INTEGRATED * PANEL. PERR IS USED DURING THE NOISE TEST. * PF1 PF2 ARE THE PREVIOUS VALUES OF F1 AND F2. PF1 AND PF2 ARE USED TO * COMPUTE EPSMIN. * PHI CONTAINS THE DIFFERENCE LINES. PHI IS DIMENSIONED AND * EQUIVALENCED WITH PHIT AFTER THE DESCRIPTION OF INTERNAL * VARIABLES. * PHIS IS A TWO-DIMENSIONAL ARRAY EQUIVALENCED TO PHI AND PHIT. IT * IS USED TO ACCESS THE DIFFERENCE LINES USING AN INDEX. */ /* PHISUM IS THE SUM OF 3 MIDDLE ENTRIES OF PHI. * PHIT IS THE BACKWARD DIFFERENCE LINE. SEE PHI. */ /* PHTSUM IS THE SUM OF 3 MIDDLE ENTRIES OF PHIT. * PX IS THE PREVIOUS VALUE OF THE OFFSET FROM THE ENDS OF THE * CURRENT INTERVAL. PX IS USED TO COMPUTE EPSMIN. * RE IS THE RATIO OF THE ERROR COMMITTED WITH THE CURRENT * QUADRATURE FORMULA TO THAT COMMITTED WITH THE PREVIOUS * QUADRATURE FORMULA. RE IS USED TO DETECT CONVERGENCE AND * TO AVOID FALSE CONVERGENCE, AND TO DETECT NOISE. * RECORR IS USED TO CORRECT RE WHEN NSUB IS NOT ZERO. ONLY THE FIRST * AND THIRD ELEMENTS OF RECORR ARE USED EXPLICITLY. */ /* RELEPS IS THE ABSOLUTE ERROR TOLERANCE, BASED UPON THE RELATIVE * RELATIVE TOLERANCE REQUESTED AND THE BEST AVAILABLE ESTIMATE * OF THE ANSWER. * RELOBT IS THE PRECISION RELATIVE TO THE OBTAINABLE PRECISION. * THE ESTIMATED ERROR IS COMPARED TO * (EPSMIN/AACUM)**RELOBT*AACUM. * RELOOK IS USED TO FORCE EXAMINATION OF THE DIFFERENCE LINES, IF * THE METHOD IS TO ACCEPT AN ANSWER BUT RE*REP .GT. RELOOK. */ /* RELTOL IS THE INDEX INTO WORK OF THE PARAMETERS FOR ABSOLUTE, LOCAL * RELATIVE, AND GLOBAL RELATIVE TOLERANCE PARAMETERS. * REPROD IS THE PRODUCT OF ALL VALUES OF RE FOR THE CURRENT PANEL. * REPROD IS NOT ALLOWED TO BE GREATER THAN 1.0. * REP IS THE PREVIOUS VALUE OF RE. IT IS USED SIMILARLY TO RE. * RESULT IS THE CUMULATIVE INTEGRAL FOR THE CURRENT MAJOR SUBDIVISION. * REVERS IS THE INDEX OF THE ENTRY IN THE OPTION VECTOR (IOPT) TO * USE TO CONTROL REVERSE COMMUNICATION. * RNDC IS THE NOISE OR ROUND-OFF LEVEL, INITIALLY EQUAL TO FER. * RNDC IS CHANGED IF NOISE IS DETECTED. * ROUNDF IS SET TO INDICATE THAT ROUND-OFF IS LIMITING THE ERROR, OR * THAT THE INTEGRAND IS BEING TREATED AS PURE NOISE WITH RESPECT * TO THE REQUESTED ERROR TOLERANCE. * S IS A SCRATCH VARIABLE. * SEARCH IS USED TO INDICATE WHETHER POINTS NOT USED BY THE INTEGRATION * FORMULAE ARE BEING ADDED TO THE DIFFERENCE LINES. * 1 - POINTS NOT BEING ADDED * 2 - SEARCHING FOR INITIAL INTERVAL SIZE * 3 - MAKING SURE AN END-POINT JUMP IS NOT INSIDE THE INTERVAL * 4 - SEARCHING FOR A DISCONTINUITY * 5 - CHANGING FROM SEARCH=2 OR 3 TO SEARCH=4 * 6 - NSUB WAS JUST REDUCED BECAUSE WE SAW JUMPS. * NSUB CANNOT BE INCREASED WITHOUT DOING A TYPE-2 SEARCH. * NEGATIVE - INITIALIZE THE SEARCH (TRANSIENT STATE). * SMIN IS THE MINIMUM STEPSIZE TO USE FOR A PANEL BEGINNING AT SX. */ /* SPACE RESERVES SOME SPACE IN /SINTC/ IN CASE IT IS NEEDED. THIS * IS DONE SO USERS' CODE WON'T NEED TO CHANGE. * START STORES THE 'LEFT' BOUNDARY OF A PART. * STEP STORES THE SIZE OF THE FIRST STEP TO USE ON A MAJOR * SUBDIVISION. THE SIGN OF STEP IS THE SIGN TO USE FOR * ACCUMULATION OF THE ANSWER INTO RESULT(PART). * SUM IS THE CENTER OF THE CURRENT INTERVAL. * T IS THE TRANSFORMED INDEPENDENT VARIABLE. SEE NSUB. * TA TB ARE THE INTERVAL ON WHICH A VARIABLE TRANSFORMATION * IS DEFINED. SEE NSUB. THE SIGN OF TB IS ALSO USED TO * INDICATE THE DIRECTION OF STEP ACCUMULATION. * TALOC IF NON-ZERO IS THE INDEX IN WORK OF THE LOCATION OF A * SINGULARITY. * TASAVE IS THE VALUE TO USE FOR TA WHEN RESTARTING PART 1 * AFTER SUBDIVIDING. * TDECR IS AN ARITHMETIC STATEMENT FUNCTION DEFINED BELOW. */ /* TEND IS THE POINT AT WHICH USAGE OF A T**N SUBSTITUTION SHOULD BE * DISCONTINUED. * THROW IS USED DURING THE SEARCH FOR A DISCONTINUITY TO INDICATE * WHICH END OF THE DIFFERENCE LINES SHOULD BE DISCARDED. */ /* TINCR IS AN ARITHMETIC STATEMENT FUNCTION DEFINED BELOW. */ /* TLEN IS THE TOTAL LENGTH OF THE UN-INTEGRATED PART. TLEN IS USED * TO PARTION ERROR ALLOTTMENTS. * TP TP1 ARE TEMPORARY VARIABLES. */ /* TPS IS A TEMPORARY VARIABLE. * WHERE IS USED AS A COMPUTED GO TO INDEX. * WHERE2 IS USED AS A COMPUTED GO TO INDEX. * WORRY IS USED TO STORE THE CAUTION POINTS FOR EACH PART. * X DURING THE QUADRATURE STEP, X IS THE DISTANCE FROM THE ENDS OF * THE INTERVAL, AT WHICH ABSCISSAE THE INTEGRAND IS EVALUATED. * X IS ALSO USED ELSEWHERE FOR TEMPORARY STORAGE. * XEPS IS THE MINIMUM PART OF THE REMAINING ERROR TOLERANCE TO * REQUIRE ON THE CURRENT INTERVAL. */ /* XJ IS THE MAGNITUDE OF THE LARGEST JUMP IN THE BACKWARD * DIFFERENCE LINE. * XJP IS THE MAGNITUDE OF THE LARGEST JUMP IN THE FORWARD * DIFFERENCE LINE. * XJUMP IS A THRESHOLD FOR JUMPS IN THE DIFFERENCE LINE. * XJUMPS IS THE VALUE OF XJUMP TO USE DURING A SEARCH. */ /* XSTEP IS USED TO INCREASE THE CURRENT EFFECTIVE STEPSIZE AFTER * SUCCESSFULLY INTEGRATING AN INTERVAL. THE CURRENT EFFECTIVE * STEPSIZE IS USED TO COMPUTE THE NEXT STEPSIZE. XSTEP IS * READ-ONLY, INDEXED BY MAX(AIMFOR-K+3,1). */ /* XT THE ABSCISSAE OF THE QUADRATURE FORMULA. XT ARE USED TO * COMPUTE THE DIFFERENCE LINE. * X1 X2 ARE TEMPORARY STORAGE. * Z ARE ADJUSTABLE PARAMETERS USED FOR DEVELOPMENT. */ /* ZL1 IS THE ARGUMENT OF ARITHMETIC STATEMENT FUNCTIONS. */ /* ***** COMMON STORAGE ****************************************** * * COMMON /SINTNC/ CONTAINS VARIABLES NOT SEPARATELY SAVED FOR * EACH DIMENSION OF A MULTIPLE QUADRATURE. COMMON /SINTC/ * CONTAINS VARIABLES THAT MUST BE SAVED FOR EACH DIMENSION OF THE * QUADRATURE. THE VARIABLES IN EACH COMMON BLOCK ARE STORED IN THE * ORDER - ALWAYS DOUBLE, DOUBLE IF DOUBLE PRECISION PROGRAM, DOUBLE * IF DOUBLE PRECISION PROGRAM AND EXPONENT RANGE OF DOUBLE AND * SINGLE VERY DIFFERENT, SINGLE, INTEGER, LOGICAL. A PAD OF LOGICAL * VARIABLES IS INCLUDED AT THE END OF /SINTC/. THE DIMENSION OF * THE PAD MAY NEED TO BE VARIED SO THAT NO VARIABLES BEYOND THE END * OF THE COMMON BLOCK ARE ALTERED. * * DECLARATIONS OF COMMON /SINTNC/ VARIABLES. * */ /* DECLARATIONS OF COMMON /SINTC/ VARIABLES. * *--D Next line special: S => D, X => Q, D => D, P => D */ /* 139 $.TYPE.$ VARIABLES */ /* Note XT, FT, and PHI above are last, because they must be in adjacent * locations in SINTC. * 30 $DSTYP$ VARIABLES */ /* 29 INTEGER VARIABLES */ /* 11 TO 18 LOGICALS (7 ARE PADDING). */ /* THE COMMON BLOCKS. * */ /* 1 2 3 4 5 6 7 8 * 9 10 11 12 13 1 2 3 * 4 (2,2) 8 9 10 11 12 13 14 * 15 16 17 18 19 20 */ /* 1 2 (4) 6 7 8 9 10 11 (2) * 13 (2) 15 16 17 (2) 19 20 (24) 44 * 45 46 47 48 49 50 51 (6) * 57 (2) 59 (2) 61 62 63 64 65 * 66 (2) 68 69 70 71 72 * 73 (17) 90 (17) 107 (34) */ /* 141 142 143 144 145 146 * 147 148 149 150 (2) 152 153 * 154 (2) 156 157 158 159 160 * 161 162 163 * 164 165 166 167 168 169 */ /* 170 171 172 * 1 2 3 4 5 6 7 8 */ /* THE VARIABLES HERE DEFINE THE MACHINE ENVIRONMENT. ALL ARE SET * IN DINTOP. THE MEANING ATTACHED TO THESE VARIABLES CAN BE * FOUND BY LOOKING AT THE DEFINITIONS IN DINTOP. */ /* ***** EQUIVALENCE STATEMENTS ******************************* * */ /* ***** DATA STATEMENTS *************************************** * */ /* K-KAIMT = -4 -3 -2 -1 0 */ /* MULTIPLY DELTA BY Z(1) IN INITIAL INTERVAL SEARCH * DELTA.LE.Z(2)*ABSDIF IN INITIAL INTERVAL SEARCH MEANS DONT GO ON * DELTA.LT.Z(3)*ABSDIF IF K.LT.6 AS FOR Z(2) * DELTA.GT.Z(4)*ABSDIF MEANS KMAX=5 AFTER INITIAL INTERVAL SEARCH * * IN THE COMMENTS BELOW, F(K,I) REFERS TO THE FUNCTION VALUE * COMPUTED FOR THE I'TH NODE OF THE K'TH FORMULA. * P Points * Indices Use (C = Correction Coefficient N&W = Nodes and Weights) * 1 3 C for F(1,1) * 2- 3 N&W for F(2,1) * 4- 5 7 C for F(1,1) and F(2,1) * 6- 9 N&W for F(3, 1-2) * 10- 13 15 C for F(1,1), F(2,1), F(3,1-2) * 14- 21 N&W for F(4, 1-4) * 22- 29 31 C for F(1,1), F(2,1), F(3,1-2), F(4,1-4). * 30- 45 N&W FOR F(5,1-8). * 46- 61 63 C for F(1,1), F(2,1), F(3,1-2), F(4,1-4), F(5,1-8). * 62- 93 N&W FOR F(6,1-16). * 94-105 127 C for F(3,1), F(4,1-2), F(5,1-3), F(6,1-6). * 106-169 N&W FOR F(7,1-32). * 170-177 256 C for F(4,1), F(5,1), F(6,1-2), F(7,1-4). * 178:305 N&W FOR F(8,1-64). * * CORRECTIONS, NODES, AND WEIGHTS FOR THE 3-POINT FORMULA. * *++ Save data by elements if ~.C. */ /* ***** STATEMENT FUNCTIONS ********************************* * * TDECR IS USED TO TRANSFORM AN ABSCISSA FROM THE CURRENT COORDINATE * SYSTEM TO ONE IN WHICH NSUB IS DECREMENTED BY A FACTOR OF 2. * TDECR(ZL1)=TA+(ZL1-TA)*((ZL1-TA)/TB) */ #define TDECR(zl1) ((float)(sintc.ta*(1.0e0 + sintc.ta/sintc.tb) + \ (zl1)*(((zl1) - sintc.ta)/sintc.tb - sintc.ta/sintc.tb))) /* TINCR IS USED TO TRANSFORM AN ABSCISSA FROM THE CURRENT COORDINATE * SYSTEM TO ONE IN WHICH NSUB IS INCREMENTED BY A FACTOR OF 2. */ #define TINCR(zl1) ((float)(sintc.ta + signf( sqrtf( fabsf( sintc.tb*\ ((zl1) - sintc.ta) ) ), sintc.tb ))) /* ***** PROCEDURES ****************************************** * */ if (sintc.where != 0) goto L_2900; /* CHECK FOR TRIVIAL CASE * */ sintc.part = 1; sintc.releps = 0.0e0; L_10: Errt[1] = 0.0e0; /*--D Next line special: S => D, X => Q */ Result[1] = 0.0e0; Xcdobt[1] = FALSE; sintc.edue2a = 0.0e0; sintc.edue2b = 0.0e0; if (sintnc.ainit == sintnc.binit) goto L_2320; /* SETUP FOR COMPLETE INTERVAL. * SET INITIAL VALUES. * */ Iflag[1] = 0; sintnc.xj = 0.e0; sintnc.xjp = 0.e0; Fat[1] = 0.e0; Fat[2] = 0.e0; sintc.errc = 0.0e0; Errt[2] = 0.0e0; Start[1] = sintnc.ainit; End[1] = sintnc.binit; sintc.epsr = fmaxf( sintc.epso, fabsf( sintc.releps ) ); sintc.epsmax = xeps*sintc.epsr; sintc.nsub = 0; sintc.ta = Start[1]; sintc.tlen = fabsf( sintnc.binit - sintnc.ainit ); Step[1] = 1.0e0; Worry[1] = End[1]; sintc.init = TRUE; sintc.did1 = FALSE; sintc.havdif = FALSE; sintc.eps = sintc.epsr; if (sintc.taloc != 0) { sintc.i = labs( sintc.taloc ); sintc.ta = Work[sintc.i]; sintc.nsub = 2; if (sintc.taloc < 0) sintc.nsub = 4; if (signf( 1.0e0, sintc.ta - Start[1] )*(End[1] - sintc.ta) > 0.0e0) { if (sintc.ta != Start[1]) { /* TA IS INSIDE THE INTERVAL. SUBDIVIDE IMMEDIATELY. */ sintc.part = 2; Step[2] = -1.0e0; End[2] = Start[1]; Worry[2] = End[2]; /*--D Next line special: S => D, X => Q */ Result[2] = 0.0e0; Xcdobt[2] = FALSE; Start[1] = sintc.ta; Start[2] = sintc.ta; sintc.tb = End[2] - sintc.ta; Step[1] = sintsm( sintc.ta ); sintc.nsubsv = sintc.nsub; sintc.tasave = sintc.ta; sintc.eps = sintc.epsr*fabsf( End[2] - Start[2] )/sintc.tlen; } } else { /* TA IS OUTSIDE THE INTERVAL. */ if (fabsf( sintc.ta - Start[1] ) >= fabsf( sintc.ta - End[1] )) { /* TA IS NEARER BINIT. TURN THE INTERVAL AROUND. */ sintc.x = End[1]; End[1] = Start[1]; Start[1] = sintc.x; Step[1] = -Step[1]; } if (sintc.ta != Start[1]) { sintc.tb = End[1] - sintc.ta; Start[1] = TINCR( Start[1] ); if (sintc.nsub == 4) Start[1] = TINCR( Start[1] ); } } } L_40: sintc.rndc = sintnc.fer; sintc.epnoiz = sintec.eninf; *fatas = FALSE; sintc.discf = 0; sintc.dischk = 0; sintc.fsaved = FALSE; sintc.kmax = sintnc.kmaxf; sintc.kaimt = min( max( 5, sintnc.kmaxf ), maxk - 2 ); *alocal = Start[sintc.part]; sintc.tb = End[sintc.part] - sintc.ta; /* GO COMPUTE MINIMUM STEPSIZE */ L_50: sintc.delmin = sintsm( *alocal ); L_60: *blocal = End[sintc.part]; sintc.tend = *blocal; sintc.iend = TRUE; L_70: *fatbs = FALSE; sintc.kmin = 0; goto L_80; L_75: sintc.eps = sintc.epsr*sintc.delta/sintc.tlen; sintc.init = FALSE; /* INITIALIZE VALUES NECESSARY TO INTEGRATE ONE PANEL. * */ L_80: Local[3] = *alocal; Local[4] = *blocal; if (sintc.nsub != 0) { Local[3] = TDECR( Local[3] ); Local[4] = TDECR( Local[4] ); if (sintc.nsub != 2) { Local[3] = TDECR( Local[3] ); Local[4] = TDECR( Local[4] ); } } sintc.delta = *blocal - *alocal; sintc.diff = 0.5e0*sintc.delta; sintc.absdif = fabsf( sintc.diff ); sintc.px = sintc.diff; sintc.delta = fabsf( sintc.delta ); sintc.sum = *alocal + 0.5e0*(*blocal - *alocal); sintc.epss = fmaxf( sintc.eps, sintc.epsmax ); sintc.did1 = Start[sintc.part] != *alocal; sintc.roundf = FALSE; sintc.re = 0.2e0; sintc.reprod = 1.0e0; sintc.ip = 1; sintc.iold = 0; sintc.k = 2; if (sintnc.iprint > 1) sinto( 1, work ); /* APPLY 1-POINT GAUSS FORMULA (MIDPOINT RULE). * */ sintc.abscis = sintc.sum; sintc.where = 1; goto L_2850; L_90: Funct[1] = sintnc.fncval; sintc.pf1 = Funct[1]; sintc.pf2 = Funct[1]; sintc.pacum = Funct[1]*sintc.delta; sintc.paacum = fabsf( sintc.pacum ); sintc.erri = sintc.pacum; sintc.err = sintc.erri; sintc.pepsmn = fabsf( sintc.pacum*sintc.rndc ); if (sintnc.fea != 0) sintc.pepsmn += fabsf( sintc.delta*sintc.t*Work[sintnc.fea] ); /* GO ON TO NEXT FORMULA. * */ sintc.acum = Funct[1]*P[1]; sintc.aacum = fabsf( sintc.acum ); sintc.epsmin = 0.0e0; goto L_130; L_95: sintc.dischk = max( sintc.dischk, 1 ); L_100: sintc.k += 1; sintc.pacum = sintc.acum; sintc.pepsmn = sintc.epsmin; sintc.epsmin = sintec.enzer; /*--D Next line special: S => D, X => Q */ sintc.acum = 0.0e0; sintc.aacum = 0.0e0; /* COMPUTE CONTRIBUTION FROM PREVIOUSLY CALCULATED FUNCTION VALUES. * */ sintnc.kk = Korect[sintc.k - 2]; L_110: sintc.j = sintnc.kk%400; sintnc.kk /= 400; sintc.j1 = sintc.j%20; sintc.j2 = sintc.j/20; for (sintc.j = sintc.j1; sintc.j <= sintc.j2; sintc.j++) { sintc.ip += 1; sintc.acum += P[sintc.ip]*Funct[sintc.j]; } if (sintnc.kk != 0) goto L_110; /* COMPUTE CONTRIBUTION FROM NEW FUNCTION VALUES. * */ L_130: sintc.inew = sintc.iold + 1; sintc.iold += sintc.inew; sintc.j2 = Fstore[sintc.k - 1]%100; sintc.j1 = Fstore[sintc.k - 1]/100; sintc.j = sintc.inew; L_140: sintc.x = P[sintc.ip + 1]*sintc.diff; sintc.x1 = *blocal - sintc.x; sintc.x2 = *alocal + sintc.x; sintc.abscis = sintc.x1; sintc.errf = 0.0e0; sintc.where = 2; if (sintc.j == sintc.inew) goto L_2840; /* THE FIRST FUNCTION VALUE REQUESTED FOR EACH FORMULA IS THE * ONE NEAREST THE EDGE OF THE PANEL. */ goto L_2850; L_150: sintc.f1 = sintnc.fncval; if (sintnc.fea != 0) sintc.errf = fabsf( sintc.t*Work[sintnc.fea] ); sintc.abscis = sintc.x2; sintc.where = 3; if (sintc.j == sintc.inew) goto L_2840; goto L_2850; L_160: sintc.f2 = sintnc.fncval; if (sintnc.fea != 0) sintc.errf += fabsf( sintc.t*Work[sintnc.fea] ); /* UPDATE TOLERANCE IF A FUNCTION VALUE WAS REQUESTED AT THE EDGE OF * THE INITIAL INTERVAL. */ if (sintc.j == sintc.inew) { neweps = FALSE; if (*alocal == sintnc.ainit) { sintc.edue2a = fabsf( sintc.f1*sintc.errina ); neweps = TRUE; } else if (*alocal == sintnc.binit) { sintc.edue2b = fabsf( sintc.f1*sintc.errinb ); neweps = TRUE; } if (*blocal == sintnc.ainit) { sintc.edue2a = fabsf( sintc.f2*sintc.errina ); neweps = TRUE; } else if (*blocal == sintnc.binit) { sintc.edue2b = fabsf( sintc.f2*sintc.errinb ); neweps = TRUE; } if (neweps) { sintc.eps = sintc.edue2a + sintc.edue2b; sintc.eps = fmaxf( sintc.epsr - sintc.eps, 0.1e0*(Errt[1] + Errt[2] + sintc.eps) )*sintc.delta/sintc.tlen; } } /* COMPUTE EPSMIN */ sintc.x1 = (fabsf( sintc.x1*(sintc.f1 - sintc.pf1) ) + fabsf( sintc.x2* (sintc.f2 - sintc.pf2) ))*sintec.emeps; sintc.x1 /= fabsf( sintc.x - sintc.px ); sintc.px = sintc.x; sintc.pf1 = sintc.f1; sintc.pf2 = sintc.f2; sintc.ip += 2; sintnc.s = fabsf( sintc.f2 ); if (sintc.nsub != 0) sintnc.s = fmaxf( sintnc.s, fabsf( *answer*(sintc.ta/(sintc.x2 - sintc.ta)) ) ); sintc.epsmin += P[sintc.ip]*(sintc.x1 + (fabsf( sintc.f1 ) + sintnc.s)* sintc.rndc + sintc.errf); sintc.aacum += P[sintc.ip]*(fabsf( sintc.f1 ) + fabsf( sintc.f2 )); sintc.f2 += sintc.f1; sintc.acum += P[sintc.ip]*sintc.f2; if (sintc.j1 <= sintc.j2) { Funct[sintc.j1] = sintc.f2; sintc.j1 += 1; } if (sintc.j <= 7) Funct[17 + sintc.j] = sintc.f1; sintc.j += 1; if (sintc.j <= sintc.iold) goto L_140; /*--D Next line special: S => D, X => Q */ sintc.acum = sintc.absdif*sintc.acum + 0.5e0*sintc.pacum; sintc.aacum = fmaxf( sintec.edelm1, sintc.absdif*sintc.aacum + 0.5e0*sintc.paacum ); sintc.epsmin = fmaxf( sintec.edelm1, sintc.absdif*sintc.epsmin + 0.5e0*sintc.pepsmn ); sintc.epsr = fmaxf( sintc.epsr, sintc.epsmin ); /* CHECK FOR CONVERGENCE. * */ sintc.ep = fmaxf( fabsf( sintc.erri ), 0.125e0*sintc.epsmin ); sintc.erri = sintc.acum - sintc.pacum; sintc.err = fabsf( sintc.erri ); sintc.extra = 0.0e0; sintc.rep = sintc.re; if (sintc.init) { if (sintnc.reltol != 0.0e0) { if (sintc.releps <= 0.0e0) { sintc.releps = fabsf( sintc.acum ) - sintc.err; sintc.releps = -fabsf( Work[sintnc.reltol + 2] )*fmaxf( sintc.releps, 0.0e0 ); sintc.epsr = fmaxf( sintc.eps, fabsf( sintc.releps ) ); sintc.epsmax = xeps*sintc.epsr; sintc.epss = fmaxf( sintc.epsmax, sintc.epsr ); } } } sintc.eps = fmaxf( sintc.epss, powf(sintc.epsmin/sintc.aacum,sintnc.relobt)* sintc.aacum ); sintc.re = sintc.err/sintc.ep; sintnc.tps = sintc.re + sintc.re; if (sintc.nsub != 0) { sintnc.tps = 0.25e0 - 1000.0e0*powif(fabsf( (sintc.ta - Local[3])/ (Local[4] - Local[3]) ),sintc.nsub); sintnc.tps = 2.0e0*fmaxf( sintc.re, Recorr[sintc.nsub - 1]* sintnc.tps ); } sintc.reprod = fminf( sintc.re*sintc.reprod, 1.0e0 ); sintc.fail = TRUE; sintc.search = 1; sintc.where = 2; sintc.havdif = sintc.havdif && sintc.k > 4; sintc.j = min( 0, sintc.k - sintc.kaimt ); if (sintc.j < 0) sintnc.tps = fmaxf( sintnc.tps, 2.0e0*fminf( 5.0e0*sintnc.tps, fminf( sintc.rep, 0.125e0 ) ) ); sintnc.tp = sintc.epsmin*(Fudge[sintc.j + 5]/sintc.rndc); if (sintc.err > sintnc.tp) { /* IF ERROR IS LARGE RELATIVE TO THE OBTAINABLE PRECISION, INCREASE * THE ERROR ESTIMATE. (THIS TEST IS CONSERVATIVE WHEN K IS SMALL) */ tp1 = 0.25e0/fmaxf( 1.0e0, 25.0e0*sintc.rep*sintc.rep ); L_190: sintc.err *= 4.0e0; sintnc.tp *= 32.0e0; sintnc.tps = sintnc.tps*(1.0e0 + sintnc.tps)/(tp1 + sintnc.tps); if (sintc.err > sintnc.tp) goto L_190; /* END OF INCREASING ERROR ESTIMATE WHEN RELATIVE ERROR IS LARGE. */ if ((sintc.err*sintc.rndc) >= sintc.epsmin) { if (sintc.err > sintc.eps) switch (IARITHIF(sintc.dischk)) { case -1: goto L_220; case 0: goto L_220; case 1: goto L_230; } /* GOT ACCURACY REQUESTED, INTEGRAND ASSUMED TO BE ALL NOISE. * DO NOT BLINDLY ACCEPT RESULT IF ON INITIAL INTERVAL. */ switch (IARITHIF(sintc.k - 2)) { case -1: goto L_220; case 0: goto L_220; case 1: goto L_230; } } } sintc.err *= sintnc.tps; /* DO NOT ACCEPT RESULT WITHOUT LOOKING AT DIFFERENCE LINE IF * WORKING ON THE INITIAL INTERVAL OR THE LEFT BOUNDARY WAS A * STOP POINT OR AN END-POINT SINGULARITY, OR IF RE*REP .GT. RELOOK. */ if (sintc.dischk > 0) goto L_230; if (sintc.err > sintc.eps) goto L_220; if (sintc.k != sintc.kaimt) goto L_230; if (sintc.dischk < 0) goto L_230; if (sintc.re*sintc.rep > relook) goto L_230; if (sintc.init) goto L_230; if (sintc.re >= 1.0e0) goto L_220; sintc.where = 3; goto L_230; L_220: sintc.where = 1; L_230: if (sintnc.iprint > 2) sinto( 2, work ); sintc.ep = sintc.errc; sintc.errc = sintc.err; sintc.erri = fabsf( sintc.erri ); switch (sintc.where) { case 1: goto L_240; case 2: goto L_250; case 3: goto L_2250; } L_240: if (sintc.k < min( sintc.kmax, sintc.kaimt )) goto L_100; /* THE ERROR IS NOT SMALL ENOUGH, OR IT SEEMS SMALL ENOUGH BUT WE * HAVE A REASON TO BE SUSPICIOUS. * */ L_250: if (sintc.havdif) goto L_500; /* COMPUTE ABSCISSAE AND UNSCRAMBLE FUNCTION VALUES INTO * USABLE ORDER. */ Ft[9] = Funct[1]; Xt[9] = sintc.sum; sintc.j = 1; sintc.l = 8; sintnc.kk = min( sintc.k, 4 ); for (sintc.j1 = 2; sintc.j1 <= sintnc.kk; sintc.j1++) { sintnc.inc = sintc.l; sintc.l /= 2; for (sintc.i = sintc.l, _do0=DOCNT(sintc.i,7,_do1 = sintnc.inc); _do0 > 0; sintc.i += _do1, _do0--) { sintc.j += 1; sintc.j2 = Ih[sintc.i]; beta = sintc.diff*P[sintc.j2]; Xt[sintc.i + 1] = *alocal + beta; Xt[17 - sintc.i] = *blocal - beta; Ft[sintc.i + 1] = Funct[sintc.j] - Funct[16 + sintc.j]; Ft[17 - sintc.i] = Funct[16 + sintc.j]; } } sintc.endpts = 1; if (sintc.k != 2) goto L_300; sintc.endpts = 3; /* COMPUTE F AT (OR NEAR) THE ENDS OF THE PANEL. */ sintc.where = 4; sintc.where2 = 1; goto L_2810; L_280: Ft[1] = *fata; Xt[1] = sintc.x; sintc.where2 = 2; goto L_2810; L_290: Ft[17] = *fatb; Xt[17] = sintc.x; L_300: sintc.where = 0; /* WHERE=0 TELLS SINTDL TO FORM AND ANALYZE DIFFERENCE LINES. */ goto L_370; /* A NEW FUNCTION VALUE HAS BEEN COMPUTED. */ L_310: sintc.where = 3; /* NSUB HAS BEEN CHANGED OR THE DIFFERENCE LINES HAVE BEEN OTHERWISE * MUTILATED. UPDATE THE DIFFERENCE LINES. */ L_350: if (sintnc.kdim == 1) sintc.errf = 0.0e0; if (sintnc.fea != 0) sintc.errf = fabsf( sintc.t*Work[sintnc.fea] ); /* WHERE TELLS SINTDU WHETHER TO ADD A FUNCTION VALUE IN THE * MIDDLE, ON THE ALOCAL END OR THE BLOCAL END, OR TO RECOMPUTE THE * DIFFERENCE LINES. */ sintdu(); L_360: sintc.where = 1; /* WHERE .NE. 0 TELLS SINTDL TO ANALYZE DIFFERENCE LINES THAT * HAVE ALREADY BEEN FORMED. */ L_370: sintdl( work ); L_500: konvrg = *ic1 + sintc.lendt - *ic2; /* EXAMINE THE ISTOP ARRAY. * */ L_510: jumps = 1; if (sintnc.istop[1][1] + 18 - sintnc.istop[0][1] == 0) goto L_530; jumps = 4; if (sintnc.istop[1][1]*(18 - sintnc.istop[0][1]) != 0) goto L_530; jumps = 3; if ((sintc.lendt - sintnc.istop[0][1])*(1 - sintnc.istop[1][1]) == 0) jumps = 2; goto L_530; /* IF WE WERE NOT DOING A JUMP SEARCH AND FOUND JUMPS, START ONE. */ L_520: sintc.esold = -1.0e0; sintc.count = 0.0e0; sintc.xjump = xjumps; sintc.search = 5; /* SET A CAUTION POINT. */ Worry[sintc.part] = Xt[sintc.lendt]; if (sintc.part == 2 && Step[1]*Step[2] < 0.0e0) Worry[1] = Xt[1]; L_530: ; /* DO BLOCK */ sintc.j1 = sintnc.istop[0][1]; sintc.j2 = sintnc.istop[1][1]; sintc.err = sintc.errc; switch (sintc.search) { case 1: goto L_540; case 2: goto L_2150; case 3: goto L_770; case 4: goto L_860; case 5: goto L_800; case 6: goto L_535; } /* One hopes SEARCH==6 and JUMPS==4 can't happen */ L_535: switch (jumps) { case 1: goto L_2160; case 2: goto L_620; case 3: goto L_970; case 4: goto L_520; } L_540: sintc.count = 0.0e0; if (jumps == 1) goto L_1920; if (sintc.k <= 2) goto L_95; if (sintc.delta < sintc.delmin + sintc.delmin) goto L_2170; if (sintc.dischk == 0) sintc.dischk = -1; if (jumps == 3) goto L_1920; if (jumps != 2) goto L_790; /* 88888888888 * END POINT JUMP. * */ if ((sintc.ta - *alocal)*signf( 1.0e0, *blocal - sintc.ta ) > 0.0e0) switch (IARITHIF(sintc.endpts - 2)) { case -1: goto L_2750; case 0: goto L_2760; case 1: goto L_560; } L_560: ; if (sintc.nsub != 0) { /* NSUB CAN ONLY BE NON-ZERO WHEN .NOT.DID1 IF THE JUMP WAS EXPECTED. */ if (!sintc.did1) sintnc.xj = fabsf( sintnc.xj ); } if (sintnc.xj >= 0.0e0) { if (sintnc.xjp > 0.0e0) { if (sintc.j2 == 0) goto L_630; if (sintc.re > 1.0e0) goto L_670; sintc.err += sintc.extra; if (sintc.err <= sintc.eps) goto L_750; /*?? Do we need to ask (re .lt. rep) and compute two ways? *## S=RE*(RE/REP)*ERR *## IF (S.GT.EPS) GO TO 670 */ if (sintc.re*sintc.re*sintc.err > sintc.rep*sintc.eps) goto L_670; /* IF KAIMT IS EVER PERMITTED .GT. MAXK-2, MAXK-K MUST BE TESTED. */ switch (IARITHIF(sintc.k - sintc.kaimt - 1)) { case -1: goto L_2740; case 0: goto L_2740; case 1: goto L_670; } } } switch (IARITHIF(sintc.endpts - 2)) { case -1: goto L_2750; case 0: goto L_2760; case 1: goto L_600; } L_600: ; if (sintc.search == 1) { /* IF THE ERROR IN THE FUNCTION AT THE BOUNDARY OF THE PANEL IS * LARGER THAN THE DIFFERENCE BETWEEN THE FUNCTION AND THE * VALUE PREDICTED BY EXTRAPOLATING THE DIFFERENCE LINE, * REPLACE THE FUNCTION BY THE EXTRAPOLATED VALUE, AND DECREASE * THE ERROR BY THE CORRECTION. THIS DOESN'T AFFECT THE ERROR * IN THE INTEGRAL BECAUSE FUNCTION VALUES AT THE ENDS DON'T * CONTRIBUTE TO THE INTEGRAL. */ Endphi[2] = sintc.lendt; for (sintc.where2 = 1; sintc.where2 <= 2; sintc.where2++) { if (Fats[sintc.where2]) { if (Errat[sintc.where2] > fabsf( phis[sintc.where2 - 1][Endphi[3 - sintc.where2] - 1] )) { Fat[sintc.where2] -= phis[sintc.where2 - 1][Endphi[3 - sintc.where2] - 1]; Ft[Endphi[sintc.where2]] = Fat[sintc.where2]; Errat[sintc.where2] -= fabsf( phis[sintc.where2 - 1][Endphi[3 - sintc.where2] - 1] ); goto L_310; } } } } sintc.dischk = max( sintc.dischk, 1 ); if (sintc.re <= 1.0e0) { if (sintc.j2 == 0) goto L_630; if (sintc.err <= sintc.eps) goto L_760; } if (sintc.k < 4) goto L_100; /* END BLOCK */ L_620: if (sintc.j2 != 0) goto L_670; /* IF WE SEE A JUMP NEAR BLOCAL AND WE ARE AT THE END OF A PART, * TURN THE INTERVAL AROUND. IF WE ARE NOT AT THE END OF A * PART, TAKE HALF THE CURRENT INTERVAL AND START OVER. * */ if (sintc.search != 1) goto L_1440; L_630: ; /* DO BLOCK */ if (sintc.iend) { if (sintc.nsub == 0) goto L_640; if (fabsf( *alocal - sintc.ta ) >= 0.125e0*fabsf( sintc.tb )) { if (!sintc.did1 || sintnc.xjp <= 0.1e0) goto L_640; sintnc.istop[0][1] = 18; /* IF THE JUMP IS WEAK, PRETEND ITS NOT THERE. */ goto L_510; } } sintc.delta *= 0.5e0; goto L_2610; /* END BLOCK */ L_640: ; if (sintc.nsub != 0) { sintns( NSRB ); goto L_640; } Start[sintc.part] = *blocal; *blocal = *alocal; End[sintc.part] = *alocal; *alocal = Start[sintc.part]; sintc.ta = *alocal; sintc.tb = End[sintc.part] - sintc.ta; if (sintc.dischk <= 0) sintc.dischk = -2; sintns( NSIA ); /* TAKE HALF OF THE INTERVAL TO AVOID GETTING IN A LOOP. */ sintc.delta *= 0.5e0; sintc.tend = *alocal + signf( sintc.delta, sintc.tb ); Step[sintc.part] = -Step[sintc.part]; goto L_2610; /* IF WE SEE AN EXPECTED AND BELIEVABLE JUMP NEAR ALOCAL, WE * SHOULD INCREASE NSUB. OTHERWISE, GO INTO THE INTERVAL SEARCH. * */ L_670: ; if ((sintnc.xj > 0.0e0 && sintc.nsub < nsubmx) && (sintc.init || sintc.lendt - *ic2 >= 2)) { if ((sintc.re >= 1.0e0 && sintc.did1) && sintc.dischk != -2) { if (sintc.k < sintc.kaimt) goto L_100; } if ((sintc.ta - *alocal)*(*blocal - sintc.ta) > 0.0e0) { sintc.ta = *alocal; sintc.tb = End[sintc.part] - sintc.ta; } sintc.kaimt = min( max( 5, sintnc.kmaxf ), maxk - 2 ); sintc.kmax = sintnc.kmaxf; if (sintc.search != 1) { *blocal = Xt[sintc.lendt]; sintc.iend = FALSE; } sintns( NSIB ); if (!sintc.iend && (sintc.part != 2 || sintc.did1)) { /* SET A CAUTION POINT. */ if (fabsf( *blocal - Start[sintc.part] ) < fabsf( sintc.tend - Start[sintc.part] )) sintc.tend = *blocal; } if (sintc.search == 1) { if (*ic2 < sintc.lendt - 2 || sintc.k > 3) { sintc.delmin = sintsm( *alocal ); sintc.eps = sintc.epss; goto L_70; } *blocal = *alocal + 0.5e0*(*blocal - *alocal); sintc.delta *= 0.5e0; } sintc.iend = FALSE; sintc.delmin = sintsm( *alocal ); goto L_75; } /* WE CAN NOT INCREASE NSUB. START A SEARCH FOR THE CORRECT * INTERVAL LENGTH. */ switch (IARITHIF(sintc.endpts - 2)) { case -1: goto L_2750; case 0: goto L_2760; case 1: goto L_2155; } /* ADD POINTS AT THE END NEAR THE JUMP UNTIL WE ARE SURE THE * JUMP IS NOT INSIDE THE INTERVAL. THIS PROCESS CONTINUES UNTIL * THE JUMP VANISHES, MOVES INSIDE THE INTERVAL, OR THE ERROR * ESTIMATED ACROSS THE JUMP IS SUFFICIENTLY SMALL. IF ALOCAL = TA * (WHICH MEANS WE ARE FAIRLY CERTAIN THE SINGULARITY IS NOT IN THE * INTERVAL), NSUB IS NOT ZERO, AND THE JUMP IS EXPECTED, ACCEPT THE * ANSWER WITHOUT DOING A SEARCH. * */ L_750: switch (IARITHIF(sintc.endpts - 2)) { case -1: goto L_2750; case 0: goto L_2760; case 1: goto L_760; } L_760: sintc.search = -3; L_770: switch (jumps) { case 1: goto L_1190; case 2: goto L_780; case 3: goto L_980; case 4: goto L_870; } L_780: if (sintc.nsub == 0) goto L_980; if (sintnc.xj <= 0.0e0) goto L_980; switch (SARITHIF(*alocal - sintc.ta)) { case -1: goto L_980; case 0: goto L_2240; case 1: goto L_980; } /* JUMPS WERE NOTICED IN THE ORIGINAL DIFFERENCE LINE ON THIS PANEL. * */ L_790: if (sintc.k + jumps < 6) goto L_100; if (max( sintc.j1, sintc.j2 ) <= sintc.nsub + 2) goto L_2130; sintc.epss = sintc.eps; if (sintc.nsub != 0) goto L_810; sintc.search = -4; sintc.xjump = xjumps; goto L_940; L_800: if (sintc.nsub == 0) goto L_850; /* IF WE HAVE NOT DONE ANYTHING, REDUCE DELTA INSTEAD OF SEARCHING. * THERE WAS A REASON FOR SETTING NSUB NON-ZERO. */ L_810: if (sintc.j2 < 2) sintc.j2 = sintc.j1; sintc.j = min( sintc.j1, sintc.j2 ); if (fabsf( 2.0e0*(Xt[1] - sintc.ta) ) <= fabsf( Xt[sintc.j] - sintc.ta )) { Worry[sintc.part] = Xt[sintc.j]; sintc.delta = 0.75e0*fabsf( Xt[sintc.j] - *alocal ); sintc.dischk = 2; goto L_2610; } /* REDUCE NSUB TO ZERO BEFORE CONTINUING THE SEARCH. */ sintc.x1 = (float)( sintc.nsub ); sintc.j = 1; for (sintc.i = 1; sintc.i <= sintc.lendt; sintc.i++) { sintc.t = (Xt[sintc.i] - sintc.ta)/sintc.tb; if (sintc.nsub != 2) sintc.t = powif(sintc.t,sintc.nsub - 1); Ft[sintc.i] /= sintc.x1*sintc.t; Xt[sintc.i] = sintc.ta + sintc.t*(Xt[sintc.i] - sintc.ta); if (sintc.i != 1) { if (Xt[sintc.i] != Xt[sintc.i - 1]) { sintc.j += 1; Xt[sintc.j] = Xt[sintc.i]; Ft[sintc.j] = Ft[sintc.i]; } } } sintc.lendt = sintc.j; sintc.search = 6; /* FINISH DECREASING NSUB */ L_840: if (sintc.nsub != 0) { sintns( NSRB ); goto L_840; } goto L_310; /* ARE JUMPS STILL WITHIN RANGE OF PREVIOUS JUMPS? */ L_850: sintc.search = 4; L_860: if (jumps < 4) goto L_870; if (sintc.j1 < min( sintc.j1old, sintc.j2old )) goto L_870; if (sintc.j2 <= max( sintc.j1old, sintc.j2old )) goto L_940; L_870: if (sintc.nfeval <= sintc.nfjump + 3) goto L_2690; goto L_310; /* THE DIFFERENCE LINE HAS BEEN BALANCED. CONTINUE THE JUMP SEARCH * OR COMPUTE THE INITIAL INTERVAL. * */ L_910: ; if (sintc.search != 4) { switch (IARITHIF(jumps - 3)) { case -1: goto L_920; case 0: goto L_930; case 1: goto L_915; } L_915: switch (IARITHIF(sintc.j2 - (sintc.nsub + 2))) { case -1: goto L_930; case 0: goto L_520; case 1: goto L_520; } L_920: sintc.i = 6 - sintc.nsub/2; sintc.j = sintc.i; if (*ic1 > sintc.j) sintc.j = min( sintc.i + 1, 6 ); sintnc.kk = 4; sintc.kmax = max( sintc.kmax, 5 ); sintc.delta = fabsf( Xt[sintc.j + 1] - *alocal ); if (sintc.nsub == 0) sintc.delta = fminf( sintc.delta, fabsf( Xt[6] - *alocal ) + 2.0e0*fabsf( Xt[6] - Xt[5] ) ); sintc.delta *= Z[1]; sintnc.tps = sintec.emeps*sintc.epsmin/sintc.rndc; /* TPS DOESN'T DEPEND ON TOLERANCE; INIT. INTERVAL DOESN'T EITHER. */ if (sintc.nsub != 0) sintnc.tps = fminf( sintc.epsmin, sintnc.tps ); alpha = (sintc.delta*Errcf[2]*(fabsf( Phi[7] ) + fabsf( Phi[8] )))/ sintnc.tps; if (alpha < 1.0e0) sintnc.kk = 5; goto L_2480; L_930: ; } if (sintc.nfeval > sintc.nfjump + 1) goto L_310; if (jumps == 1) goto L_1190; /* JUMPS HAVE NOT ESCAPED. * */ L_940: if (sintc.j1 >= sintc.j2) goto L_970; L_945: ; if (sintc.j2 - sintc.j1 <= 1) { sintc.i = sintc.j1; sintc.j1 = sintc.j2; sintc.j2 = sintc.i; } else { /* JUMPS DO NOT OVERLAP, CHOOSE ONE ON WHICH TO CONVERGE. */ if (sintc.part != 1 && Step[1]*Step[2] <= 0.0e0) { sintc.j1 = 18; } else { sintc.j2 = 0; } } /* JUMPS OVERLAP NOW. */ L_970: sintc.j1old = sintc.j1; if (sintc.j1old == sintc.lendt) sintc.j1old = sintc.lendt - 1; L_980: sintc.j2old = sintc.j2; if (sintc.j2old == 1) sintc.j2old = 2; /* WE ONLY GET HERE IF WE HAVE AT LEAST 1 JUMP. */ if (sintc.j1 == 18) sintc.j1 = sintc.j2 + 1; if (sintc.j2 == 0) sintc.j2 = sintc.j1 - 1; if (min( sintnc.istop[0][1], sintc.lendt + 1 - sintnc.istop[1][1] ) < (sintc.lendt + 2)/3) goto L_1350; if (labs( sintc.search ) != 4) goto L_990; /* IF (JUMPS.NE.4) GO TO 1350 * The above test prevented getting out of the search early on * rapidly decaying integrands. */ if (sintnc.istop[0][1] < sintnc.istop[1][1]) goto L_1350; /* ESTIMATE ERROR ACROSS JUMP. * */ L_990: sintc.where = 0; L_1000: ; if (labs( sintc.search ) != 4 && labs( sintc.search ) != 3) { if (jumps >= 3) { /* IF (WHERE.EQ.0) GO TO 1350 * The above test prevented estimating error if JUMPS == 3 */ sintnc.inc = sintc.lendt; sintnc.inc2 = 1; } } else { if (jumps == 3) { sintnc.inc = min( sintc.lendt, sintc.j1old + 1 ); sintnc.inc2 = max( 1, sintc.j2old - 1 ); } else { sintnc.inc = min( sintc.lendt, sintc.j1 + 1 ); sintnc.inc2 = max( 1, sintc.j2 - 1 ); } } sintnc.tps = 0.0e0; sintc.f1 = -1.0e0; for (sintc.i = sintnc.inc2; sintc.i <= (sintnc.inc - 1); sintc.i++) { sintc.f2 = fmaxf( fabsf( Xt[sintc.i + 1] - Xt[sintc.i] ), fabsf( sintec.eepsm8* Xt[sintc.i] ) ); sintc.f2 = fabsf( (Ft[sintc.i + 1] - Ft[sintc.i])*sintc.f2 ); /* F1 IS THE ERROR COMMITTED ON THE TRAPEZOID WITH THE LARGEST ERR */ sintc.f1 = fmaxf( sintc.f1, sintc.f2 ); sintnc.tps += sintc.f2; } sintnc.tps *= 0.5e0; if (fmaxf( sintnc.xj, sintnc.xjp ) > 0.05e0) sintnc.tps *= 4.0e0; sintc.x = fabsf( 0.5e0*(Xt[sintnc.inc] - Xt[sintnc.inc2]) ); if (sintnc.inc2 >= sintnc.inc - 2) { if (sintnc.inc2 > 1) sintnc.tps += fabsf( sintc.x*(Ft[sintnc.inc2] - Ft[sintnc.inc2 - 1]) ); if (sintnc.inc < sintc.lendt) sintnc.tps += fabsf( sintc.x*(Ft[sintnc.inc + 1] - Ft[sintnc.inc]) ); } /* COMPUTE THE ALLOWABLE ERROR. * DO BLOCK */ if (sintc.where == 0) { /* THE ABSCISSAE DID NOT COALESCE DURING INSERTION OF A POINT * IN THE DIFFERENCE LINES */ sintnc.s = fmaxf( 0.25e0*sintc.eps, sintc.epsmax ); /* TRY TO GET AT LEAST 3 DIGITS IF THE JUMPS LOOK REAL: */ if (jumps != 3) sintnc.s = fminf( sintnc.s, 1.0e-3*(sintc.epsmin/sintc.rndc) ); if (sintc.search != 6) { switch (IARITHIF(labs( sintc.search ) - 3)) { case -1: goto L_1060; case 0: goto L_1080; case 1: goto L_1100; } L_1060: ; } /* ABS(SEARCH)=2 OR 6 HERE */ sintnc.tps *= 10.0e0; if (sintnc.iprint > 3) sinto( 6, work ); if (sintnc.tps <= 0.01e0*sintnc.s) switch (SARITHIF(sintc.err - sintnc.s)) { case -1: goto L_2240; case 0: goto L_2240; case 1: goto L_1150; } if (jumps >= 3) goto L_1350; if (sintc.nsub != 0) switch (SARITHIF(fabsf( Xt[4] - *alocal ) - sintc.delmin)) { case -1: goto L_1150; case 0: goto L_1150; case 1: goto L_1440; } if (fabsf( Xt[9] - *alocal ) >= 8.0e0*sintc.delmin) goto L_1440; /* SET STEPSIZE TO DELMIN. */ sintc.delta = 0.0e0; goto L_2620; /* ABS(SEARCH)=3 HERE */ L_1080: ; if (sintnc.tps <= sintc.eps) { if (sintnc.xj > 0.0e0) { if (sintnc.xjp > 0.0e0) { if (sintc.re <= 1.0e0) { if (sintc.nsub != 0) { /* AUGMENT EPSMIN DUE TO NOISE INTRODUCED BY * T**2N TRANSFORMATION. */ sintc.x1 = *alocal - sintc.ta; if (sintc.x1 == 0.0e0) sintc.x1 = Xt[2] - sintc.ta; sintc.x1 = fabsf( sintc.tb/sintc.x1 ); if (sintc.nsub != 2) sintc.x1 = powif(sintc.x1,sintc.nsub - 1); sintc.epsmin += fabsf( sintc.rndc*Ft[1]* sintc.ta*sintc.x1 ); } goto L_2240; } } } if (sintnc.istop[1][1] != 0 && sintnc.tps <= sintnc.s) goto L_1150; } L_1100: ; if (sintc.search <= 0) { if (sintc.search == -3) goto L_1120; /* THE JUMP SEARCH IS JUST STARTING. IF THE SUM OF THE * ERROR ESTIMATED BY THE QUADRATURE STEP, AND THE ERROR * ACROSS THE JUMP IS LESS THAN THE TOLERANCE FOR THE PANEL, * ACCEPT THE ANSWER. */ sintc.err = sintnc.tps + sintc.errc; if (sintc.err > sintc.eps) goto L_1120; sintc.discf = 1; sintnc.s = sintc.eps; goto L_1180; } } if (sintc.esold > 0.0e0) sintc.count += fminf( (sintc.f1 - sintc.esold)/sintc.esold, 1.0e0 ); /* END BLOCK */ L_1120: ; sintc.esold = sintc.f1; if (sintnc.iprint > 3) sinto( 6, work ); sintc.search = labs( sintc.search ); if (sintc.where == 0) { /* THE ABSCISSAE DID NOT COALESCE DURING INSERTION OF A POINT IN * THE DIFFERENCE LINES. */ if (sintc.search == 3) goto L_1350; switch (SARITHIF(sintnc.tps - sintnc.s)) { case -1: goto L_2670; case 0: goto L_2670; case 1: goto L_1350; } } if (sintc.search == 3) switch (SARITHIF(sintc.count)) { case -1: goto L_2237; case 0: goto L_2237; case 1: goto L_2190; } if (sintc.count <= 0.0e0) goto L_1150; Local[3] = Xt[sintnc.inc]; goto L_2190; /* THE ESTIMATED ERROR ACROSS THE JUMP IS SMALL ENOUGH TO ACCEPT. * ESTIMATE THE INTEGRAL USING THE COMPOSITE TRAPEZOID RULE. * THIS IS SATISFACTORY SINCE THE PRESENCE OF THE JUMP NEGATES * ANY POSSIBILITY OF HIGH-ORDER CONVERGENCE. * */ L_1150: sintc.eps = sintnc.s; L_1160: sintc.err = sintnc.tps; sintc.erri = sintc.err; /*--D Next line special: S => D, X => Q */ sintc.acum = 0.0e0; sintc.x1 = 0.0e0; sintc.x2 = 0.0e0; sintc.j = sintnc.inc - 1; /* DO BLOCK */ L_1170: sintnc.s = fabsf( Xt[sintc.j + 1] - Xt[sintc.j] ); /* COMPUTE EPSMIN ALSO */ sintc.f1 = fabsf( Ft[sintc.j] ); sintc.f2 = fabsf( Ft[sintc.j + 1] ); sintc.x1 += fabsf( (Xt[sintc.j + 1] + Xt[sintc.j])*(sintc.f2 - sintc.f1) ); sintc.x2 += sintnc.s*(sintc.f2 + sintc.f1); sintc.acum += (Ft[sintc.j + 1] + Ft[sintc.j])*sintnc.s; sintc.j -= 1; if (sintc.j >= sintnc.inc2) goto L_1170; /* END BLOCK */ sintc.epsmin = 0.5e0*(sintc.x1*sintec.emeps + sintc.x2*sintc.rndc); /*--D Next line special: S => D, X => Q */ sintc.acum *= 0.5e0; sintc.iend = FALSE; sintc.havdif = FALSE; /* SET K SMALL FOR STEPSIZE SELECTION PURPOSES. * K=2 */ L_1180: sintc.fail = FALSE; sintc.i = sintnc.inc; sintc.j = sintnc.inc2; switch (IARITHIF(sintc.search - 4)) { case -1: goto L_1240; case 0: goto L_1230; case 1: goto L_1240; } /* IF BOTH JUMPS VANISH, AND THE ORIGINAL ERROR WAS SUFFICIENTLY * SMALL, AND ENOUGH OF THE INTERVAL IS BEING EXAMINED, ACCEPT * THE ANSWER. * */ L_1190: if (sintc.errc <= sintc.eps && 10.0e0*fabsf( Xt[sintc.lendt] - Xt[1] ) > sintc.delta) goto L_2237; /* COMPUTE A NEW PANEL TO INTEGRATE. * */ sintc.dischk = 2; if (sintc.part != 1) { if (Step[1]*Step[2] <= 0.0e0) { if (sintc.j1 == 18) sintc.j1 = 0; sintc.j = max( 2, max( sintc.j2, sintc.j1 - 2 ) ); sintc.i = sintc.lendt; if (jumps != 1) sintc.i = min( sintc.lendt, sintc.j + 3 ); goto L_1240; } } if (sintc.j2 == 0) sintc.j2 = 18; sintc.i = min( sintc.lendt - 1, min( sintc.j1, sintc.j2 + 2 ) ); sintc.j = 1; if (jumps != 1) sintc.j = max( 1, sintc.i - 3 ); goto L_1250; /* GET THE BOUNDARIES OF THE SUBDIVISION. * */ L_1230: if (sintc.part == 1 || Step[1]*Step[2] > 0.0e0) goto L_1250; L_1240: sintc.l = sintc.i; sintc.i = sintc.j; sintc.j = sintc.l; sintc.x = sintnc.xj; sintnc.xj = sintnc.xjp; sintnc.xjp = sintc.x; L_1250: sintc.x1 = Xt[sintc.i]; sintc.x2 = Xt[sintc.j]; *fata = Ft[sintc.i]; *fatb = Ft[sintc.j]; *fatas = TRUE; *fatbs = TRUE; sintc.where = 5; if (!sintc.fail) { if (sintc.nsub == 0) { Discx[1] = sintc.x1; Discx[2] = sintc.x2; } else { Discx[1] = sintc.ta + (sintc.x1 - sintc.ta)*powif((sintc.x1 - sintc.ta)/sintc.tb,sintc.nsub - 1); Discx[2] = sintc.ta + (sintc.x2 - sintc.ta)*powif((sintc.x2 - sintc.ta)/sintc.tb,sintc.nsub - 1); } if (sintc.search == -4) goto L_2240; sintc.delta = fabsf( sintc.x1 - sintc.x2 ); if (sintnc.iprint > 1) sinto( 7, work ); if (sintc.search != 4) { *blocal = sintc.x2; sintc.absdif = 0.5e0*sintc.delta; if (fminf( fabsf( sintnc.xj ), fabsf( sintnc.xjp ) ) <= 1.0e-3) sintc.discf = 1; if (sintc.search == 2 && sintc.nsub != 0) sintns( NSRB ); goto L_2240; } } /* DECREASE NSUB TO ZERO. */ L_1260: if (sintc.nsub != 0) { sintns( NSRB ); goto L_1260; } if (sintc.fail) { /* TRY NOT TO CHOOSE A PANEL TO BE INTEGRATED THAT WILL LEAVE * A RIDICULOUSLY SMALL UN-INTEGRATED PIECE. * * DO BLOCK */ if (fabsf( sintc.x2 - *alocal ) >= 1.5e0*sintc.delmin) { if (fabsf( sintc.x2 - sintc.x1 ) <= 8.0e0*fabsf( sintc.x1 - *alocal )) { smin = sintsm( sintc.x1 ); if (fabsf( sintc.x1 - *blocal ) >= 1.5e0*smin) { if (fabsf( sintc.x2 - sintc.x1 ) <= 8.0e0*fabsf( sintc.x1 - *blocal )) goto L_1280; } if (sintc.part == 2) goto L_1270; if (!sintc.iend) goto L_1270; sintc.x1 = sintc.x2; } } sintc.delta = fabsf( sintc.x1 - *alocal ); *fatas = FALSE; goto L_2610; /* END BLOCK */ L_1270: ; sintc.x1 = *blocal; *fatas = FALSE; } L_1280: ; *fatbs = *fatas; /* SPLIT THE PROBLEM INTO TWO PARTS AND INTEGRATE AWAY FROM * THE MIDDLE IN BOTH DIRECTIONS. * * DO BLOCK */ if (sintc.part != 1) { sintc.i = 2; /*## EPSR=EPSR+ERRT(2) */ if (Step[1]*Step[2] <= 0.0e0) goto L_1290; } End[2] = *alocal; sintc.i = 1; /* END BLOCK */ L_1290: ; sintc.tlen = fabsf( End[2] - End[1] ); sintc.fsave = *fata; sintc.fsaved = *fatas; sintnc.s = Start[1]; Start[1] = sintc.x1; if (sintc.fail) { /* Jump disappeared. X1..X2 not yet integrated. */ Start[2] = sintc.x1; /* Why was the following ever done ??? * TA=X1+0.5e0*(X2-X1) */ sintc.ta = sintc.x1; sintc.tasave = sintc.ta; } else { /* Jump didn't disappear. X1..X2 was integrated using trapezoid */ Start[2] = sintc.x2; sintc.ta = Start[1]; sintc.tasave = Start[2]; } sintc.tb = End[2] - sintc.ta; sintc.f1 = Worry[1]; Worry[1] = Local[3 - sintc.i]; if (fabsf( sintc.f1 - sintnc.s ) >= fabsf( sintc.x1 - sintnc.s )) { if (fabsf( sintc.f1 - sintnc.s ) < fabsf( Worry[1] - sintnc.s )) Worry[1] = sintc.f1; } Worry[2] = Local[sintc.i]; sintc.delta = fabsf( sintc.x1 - sintc.x2 ); Step[2] = signf( 1.0e0, signf( 1.0e0, sintc.tb )*(sintnc.binit - sintnc.ainit) ); *alocal = sintc.x1; *blocal = sintc.x2; sintc.part = 2; sintc.did1 = FALSE; sintc.init = FALSE; sintc.nsubsv = 0; if (!sintc.fail) { if (sintnc.xj > 0.0e0) sintc.nsubsv = 2; if (sintnc.xjp > 0.0e0) sintns( NSIB ); } sintc.absdif = 0.5e0*sintc.delta; Errt[2] = 0.0e0; Result[2] = 0.0e0; Xcdobt[2] = FALSE; sintc.dischk = max( sintc.dischk, 1 ); if (!sintc.fail) { smin = sintsm( *blocal ); sintc.iend = fabsf( End[2] - *blocal ) <= smin; /* SOMEDAY, WE SHOULD USE A TRAPEZOIDAL ESTIMATE ACROSS ANY * TINY PIECE LEFT UN-INTEGRATED DUE TO THE ABOVE STATEMENT. */ if (jumps == 3) { sintc.discf = -2; } else { sintc.discf = 2; } sintnc.s = Worry[1]; if (sintnc.xj <= 0.0e0) sintnc.s = End[1]; /* SET THE RESTARTING STEPSIZE FOR PART 1. */ smin = sintsm( Start[1] ); Step[1] = signf( fmaxf( fabsf( sintnc.s - Start[1] ), smin ), Step[1] ); goto L_2250; } /* SET THE RESTARTING STEPSIZE FOR PART 1 TO THE STEPSIZE FOR PART 2. */ smin = sintsm( sintc.x1 ); Step[1] = signf( fmaxf( sintc.delta, smin ), Step[1] ); sintc.kmin = 0; goto L_2615; /* THE ESTIMATED ERROR ACROSS THE INTERVAL IS TOO LARGE. REDUCE * THE INTERVAL LENGTH AND ADD A FUNCTION VALUE IN THE NEIGHBORHOOD * OF THE JUMP. THEN RE-FORM AND RE-EXAMINE THE DIFFERENCES. * */ L_1350: sintc.i = sintc.j1; sintc.j = sintc.j2; /* DO FOREVER */ L_1360: ; /* DO BLOCK */ if (sintc.i != sintc.lendt) { if (sintc.j != 1) { if (jumps >= 3) { sintc.x1 = fabsf( Xt[sintc.i] - Xt[sintc.i + 1] ); sintc.x2 = fabsf( Xt[sintc.j] - Xt[sintc.j - 1] ); if (0.3e0*sintc.x2 >= fminf( sintc.x1, fabsf( Xt[sintc.j + 1] - Xt[sintc.j] ) )) { sintc.i = sintc.j - 1; goto L_1400; } if (0.3e0*sintc.x1 >= fminf( sintc.x2, fabsf( Xt[sintc.i] - Xt[sintc.i - 1] ) )) { sintc.j = sintc.i; sintc.i += 1; goto L_1400; } } } } switch (IARITHIF(sintc.i - sintc.j - 1)) { case -1: goto L_1390; case 0: goto L_1450; case 1: goto L_1410; } L_1390: sintc.i = sintc.j + 1; if (fabsf( Xt[sintc.j] - Xt[sintc.j + 1] ) < fabsf( Xt[sintc.j] - Xt[sintc.j - 1] )) sintc.i = sintc.j - 1; /* END BLOCK */ L_1400: ; sintc.abscis = Xt[sintc.j] + 0.25e0*(Xt[sintc.i] - Xt[sintc.j]); goto L_1490; L_1410: ; if (fabsf( sintnc.xj ) <= fabsf( sintnc.xjp )) { sintnc.xj *= 5.0e0; sintc.i -= 1; } else { sintnc.xjp *= 5.0e0; sintc.j += 1; } goto L_1360; /* END FOREVER */ L_1440: sintc.j = 1; L_1450: sintc.i = sintc.j + 1; sintc.abscis = Xt[sintc.i] + 0.5e0*(Xt[sintc.j] - Xt[sintc.i]); L_1490: sintc.l = max( sintc.i, sintc.j ); if (sintc.search <= 0) { sintc.search = -sintc.search; /* IF (SEARCH.EQ.4) THEN * DETERMINE WHETHER NODES HAVE COALESCED. WE ASSUME THAT NODES * WILL NOT COALESCE WHEN SEARCH = 2. * END IF */ } L_1500: if (sintc.abscis == Xt[sintc.i] || sintc.abscis == Xt[sintc.j]) { sintc.x1 = Xt[sintc.i] + 0.5e0*(Xt[sintc.j] - Xt[sintc.i]); if (sintc.x1 == sintc.abscis) goto L_2680; sintc.abscis = sintc.x1; goto L_1500; } L_1520: ; if (sintc.lendt > 7) { /* DECIDE HOW MUCH OF THE DIFFERENCE LINE TO THROW AWAY. * THROW INDICATES WHICH END TO TRIM - 1=ALOCAL, 2=BLOCAL, 3=BOTH. */ throw = 0; if (jumps >= 3) switch (IARITHIF(labs( sintc.search ) - 4)) { case -1: goto L_1555; case 0: goto L_1530; case 1: goto L_1555; } if (sintc.search == 2 || sintc.search == 6) switch (IARITHIF(sintc.lendt - 9)) { case -1: goto L_1720; case 0: goto L_1680; case 1: goto L_1630; } L_1530: ; if (sintnc.istop[1][1] > sintnc.istop[0][1] + 1) { /* IF JUMPS DO NOT OVERLAP, THROW AWAY A LOT. */ if (sintc.part != 1) { if (Step[1]*Step[2] <= 0.0e0) { sintc.j1 = sintnc.istop[1][1] - 4; sintc.j2 = 0; switch (IARITHIF(sintc.j1)) { case -1: goto L_1610; case 0: goto L_1610; case 1: goto L_1640; } } } sintc.j2 = sintc.lendt - sintnc.istop[0][1] - 3; sintc.j1 = 0; switch (IARITHIF(sintc.j2)) { case -1: goto L_1610; case 0: goto L_1610; case 1: goto L_1640; } } if (jumps == 1) goto L_1610; L_1555: sintc.j = sintc.j2; sintc.j2 = max( min( sintc.lendt - sintc.l, (*ic2 - sintc.j1 - 1)/2 ), 0 ); sintc.j1 = max( min( sintc.l - 1, (sintc.j - *ic1 - 1)/2 ), 0 ); if (sintc.search != 4) { if (sintc.lendt < 17) sintc.j1 = 0; } /* DO BLOCK */ switch (IARITHIF(sintc.j2 - 1)) { case -1: goto L_1580; case 0: goto L_1570; case 1: goto L_1640; } L_1570: throw = 2; L_1580: switch (IARITHIF(sintc.j1 - 1)) { case -1: goto L_1600; case 0: goto L_1590; case 1: goto L_1640; } L_1590: throw += 1; L_1600: if (throw != 0) switch (IARITHIF(throw - 2)) { case -1: goto L_1670; case 0: goto L_1680; case 1: goto L_1670; } /* END BLOCK */ L_1610: ; if (sintc.lendt < 15) goto L_1740; if (jumps != 1) { if (sintnc.istop[0][1] == 18) goto L_1680; if (sintnc.istop[1][1] == 0) goto L_1670; } switch (IARITHIF(sintc.lendt + 1 - sintc.j1old - sintc.j2old)) { case -1: goto L_1670; case 0: goto L_1660; case 1: goto L_1680; } L_1630: sintc.j1 = 0; sintc.j2 = sintc.lendt - 9; /* THROW OUT J1 NODES FROM THE HEAD AND J2 NODES FROM THE TAIL * OF THE DIFFERENCE LINES. */ L_1640: sintc.lendt += -sintc.j1 - sintc.j2; /* FORCE RE-FORMATION OF DIFFERENCES FROM XT AND FT. */ sintc.nfjump = sintc.nfeval - 10; if (sintc.j1 == 0) goto L_1740; if (sintc.search != 4) goto L_520; if (sintc.j1old != 18) sintc.j1old -= sintc.j1; sintc.j2old -= sintc.j1; sintc.l -= sintc.j1; for (sintc.i = 1; sintc.i <= sintc.lendt; sintc.i++) { Xt[sintc.i] = Xt[sintc.i + sintc.j1]; Ft[sintc.i] = Ft[sintc.i + sintc.j1]; } goto L_1740; L_1660: throw = 3; L_1670: if (sintc.search != 4) goto L_520; /* THROW OUT X(1) AND F(1), RE-FORM THE DIFFERENCE LINES. */ sintnc.kk = 1; sintnc.inc = 1; sintc.j = 1; sintc.j1 = 1; sintc.j2 = sintc.lendt; *fatas = FALSE; sintc.l -= 1; if (sintc.j1old != 0) sintc.j1old -= 1; if (sintc.j2old != 0) sintc.j2old -= 1; goto L_1690; /* THROW OUT X(LENDT) AND F(LENDT), RE-FORM THE DIFFERENCE LINES. */ L_1680: sintnc.kk = 2; sintnc.inc = -1; sintc.j = 17 + sintc.lendt; sintc.j1 = sintc.lendt; sintc.j2 = 1; *fatbs = FALSE; L_1690: beta = 1.0e0; sintc.i = sintc.j1; Phi[sintc.j] -= Phi[sintc.j + sintnc.inc]; sintc.j += sintnc.inc; sintc.i += sintnc.inc; /* DO BLOCK */ L_1700: beta *= (Xt[sintc.j1 + sintnc.inc] - Xt[sintc.i + sintnc.inc])/ (Xt[sintc.j1] - Xt[sintc.i]); Phi[sintc.j] = beta*(Phi[sintc.j] - Phi[sintc.j + sintnc.inc]); sintc.j += sintnc.inc; sintc.i += sintnc.inc; if (sintc.i != sintc.j2) goto L_1700; /* END BLOCK */ sintc.lendt -= 1; for (sintc.i = 1; sintc.i <= sintc.lendt; sintc.i++) { Phit[sintc.i] = Phit[sintc.i + 1]; if (sintnc.inc >= 0) { Xt[sintc.i] = Xt[sintc.i + 1]; Ft[sintc.i] = Ft[sintc.i + 1]; } } throw -= sintnc.kk; if (throw > 0) goto L_1680; /* THE NEXT STATEMENT OCCASSIONALLY CAUSES A DIFFERENCE LINE TO * BE SHORTENED, BUT NO NEW POINT ADDED. */ if (sintnc.kk + sintc.l == 2) goto L_360; L_1720: ; } L_1740: ; if (sintc.l > sintc.lendt) goto L_310; /* ADD A POINT IN THE INTERIOR OF THE INTERVAL. */ sintc.x = sintc.abscis; sintc.where = 6; goto L_2850; /* NO JUMPS, OR ONLY ONE INTERIOR JUMP, IN THE ORIGINAL * DIFFERENCE LINE FOR THIS PANEL. * */ L_1920: ; /* DO BLOCK */ if (sintc.re <= 1.0e0) { if ((sintc.k <= 5 && sintc.kaimt > 5) && konvrg < 4) sintc.err = fmaxf( sintc.err, sintc.erri*(sintc.erri/sintc.epsmin) ); if (sintc.err > sintc.eps) { if (sintc.k < 4) goto L_100; if (sintc.k > 5) goto L_2060; switch (IARITHIF(konvrg - 6)) { case -1: goto L_2060; case 0: goto L_2060; case 1: goto L_2130; } } /* ERROR CRITERIA HAVE APPARENTLY BEEN SATISFIED. LOOK AT THE * DIFFERENCE LINES TO AVOID BEING FOOLED BY ANSWERS FROM * SUCCESSIVE FORMULAE BEING NEARLY EQUAL BUT BOTH WRONG. * */ if ((konvrg <= sintc.lendt && sintc.k < sintc.kmin) && sintc.phisum <= 5.0e0*sintc.phtsum) goto L_100; if (!sintc.init && sintc.nsub == 0) { if (sintc.err > sintc.epsmin) { sintc.i = sintc.kaimt - 2; if (konvrg > sintc.lendt + sintc.lendt - 5) sintc.i = sintc.kaimt - 3; switch (IARITHIF(sintc.k - sintc.i)) { case -1: goto L_100; case 0: goto L_1950; case 1: goto L_1960; } L_1950: switch (IARITHIF(sintc.endpts - 2)) { case -1: goto L_2750; case 0: goto L_2760; case 1: goto L_1960; } L_1960: ; } } if (sintc.k <= 4) { sintc.extra = Errcf[sintc.k - 1]; if (sintc.k != 2) { if (konvrg < sintc.k + sintc.k - 1 + sintc.endpts/ 2) sintc.extra *= 10.0e0; } sintc.extra *= sintc.delta*(fabsf( Phi[sintc.lendt - 1] ) + fabsf( Phit[2] )); sintc.j = (sintc.k - 2)*(konvrg - 7); if (sintc.j > 20) { sintc.extra = 0.0e0; } else if (sintc.j > 0) { sintc.extra *= powif(0.1e0,sintc.j); } if (jumps >= 3) goto L_2140; if (sintc.extra >= 32.0e0*sintc.epsmin) { sintnc.s = sintc.err + sintc.extra; if (sintnc.s > sintc.eps) switch (IARITHIF(sintc.k - sintc.kaimt)) { case -1: goto L_2730; case 0: goto L_2050; case 1: goto L_2050; } sintc.err = sintnc.s; } } if (jumps == 3) goto L_750; /* DO BLOCK */ if (sintc.dischk <= 0) { if (sintc.dischk >= -1) { if (!sintc.init) goto L_2235; if (sintc.k >= 4) goto L_2235; sintc.kmax = maxk; goto L_2030; } } /* COMPUTE F AT (OR NEAR) THE ENDS OF THE PANEL, UNLESS * THERE IS A JUMP AT THE OTHER END OF THE PANEL. */ sintc.kmax = sintc.kaimt + 1; /* END BLOCK */ L_2030: ; switch (IARITHIF(sintc.endpts - 2)) { case -1: goto L_2750; case 0: goto L_2760; case 1: goto L_2235; } } /* CHECK FOR NOISE IN THE INTEGRAND. * */ L_2050: if (sintc.k < 4) goto L_100; L_2060: ; if (sintc.re >= 0.01e0) { /* DO BLOCK */ if (sintc.re < 1.0e0) { if (sintc.nsub != 0 && fabsf( (Local[3] - sintc.ta)/(Local[4] - Local[3]) ) <= 1.0e0) goto L_2130; } else if (sintc.erri < 100.0e0*sintc.epsmin) { if (4.0e0*sintc.err > sintc.epsmin) goto L_2100; if (!sintc.init && sintc.dischk <= 0) goto L_2250; switch (IARITHIF(sintc.endpts - 2)) { case -1: goto L_2750; case 0: goto L_2760; case 1: goto L_2250; } } if (sintc.reprod >= 1.0e-4) goto L_2130; if (sintc.k < sintc.kaimt && sintc.err > 100.0e0*sintc.epsmin) goto L_2130; if (sintc.erri*sintc.reprod > 0.01e0*sintc.epnoiz) goto L_2130; if (sintc.re*sintc.err < sintc.eps) goto L_2130; /* END BLOCK */ L_2100: ; switch (IARITHIF(sintc.endpts - 2)) { case -1: goto L_2750; case 0: goto L_2760; case 1: goto L_2120; } /* NOISE IN THE INTEGRAND - RESET THE ROUND-OFF CONSTANT. */ L_2120: sintc.rndc *= 4.0e0*(sintc.err/sintc.epsmin); if (sintnc.iprint > 0 && 100.0e0*sintc.err > sintc.eps) sinto( 8, work ); sintc.eps = sintc.err; sintc.roundf = TRUE; sintc.delta += sintc.delta; goto L_2250; } /* END BLOCK * NO NOISE DETECTED. */ L_2130: if (sintc.delta < sintc.delmin + sintc.delmin) goto L_2250; sintc.kmax = maxk; if (sintc.dischk > 0 || sintc.dischk < -1) sintc.kmax = sintc.kaimt + 1; /* IF WE HAVE A LOT OF CONVERGENCE, WE SHOULD MAKE IT. * IF WE DONT MAKE IT, SOMETHING FISHY IS HAPPENING. * ENDPTS IS INCLUDED BECAUSE CONVERGENCE IS EASIER WHEN THE END * POINTS HAVE BEEN ADDED TO THE DIFFERENCE LINES. */ if (konvrg >= sintc.endpts + 7 && sintc.re > 1.0e-3) sintc.kmax = 5; if (10.0e0*sintc.re*sintc.re*sintc.err < sintc.rep*sintc.eps) sintc.kmax = min( maxk, max( sintc.kmax, sintc.k + 1 ) ); if (jumps < 3) { if (sintc.k < sintc.kmax) { if (sintc.dischk >= -1) { if (fabsf( Worry[sintc.part] - Start[sintc.part] ) - sintc.delmin > fabsf( *blocal - Start[sintc.part] ) || sintc.did1) goto L_2730; if (sintc.k < sintc.kaimt) { sintnc.s = sintc.re*sintc.err; if (sintc.k < sintc.kaimt - 1) sintnc.s = sintnc.s*sintc.re/sintc.rep; if (sintnc.s <= sintc.eps) goto L_2730; } } } } /* USE THE SEARCH MECHANISM TO FIND THE PROPER INTERVAL LENGTH. * */ L_2140: if (sintc.endpts < 2) goto L_2750; sintc.endpts = 3; sintc.search = 2; L_2150: switch (jumps) { case 1: goto L_2160; case 2: goto L_620; case 3: goto L_970; case 4: goto L_2165; } L_2155: sintc.search = 2; L_2160: sintnc.inc = 4; sintnc.inc2 = 1; if (jumps != 1) goto L_990; switch (IARITHIF(*ic1 + *ic1 + sintc.nsub - 12)) { case -1: goto L_990; case 0: goto L_2690; case 1: goto L_2690; } L_2165: if (sintc.j1 < sintc.j2) goto L_945; sintc.xjump = xjumps; goto L_870; /* ABSCISSAE HAVE COALESCED. * */ L_2170: if (sintc.k < 4) goto L_100; L_2180: sintc.delta *= 100.0e0; if (sintc.re < 0.875e0) goto L_2200; if (sintc.erri < fmaxf( sintc.epsmin, sintc.pepsmn )) goto L_2200; if (sintc.err < fmaxf( sintc.epsmin, sintc.pepsmn )) goto L_2200; if (sintc.reprod <= 1.0e-4) goto L_2200; if (sintc.havdif && labs( jumps - 3 ) != 1) goto L_2200; L_2190: ; sinto( 9, work ); Iflag[1] = 5 + sintnc.kdim; Work[1] = sintec.eminf; *answer = *alocal; goto L_2340; L_2200: ; if (sintc.where < 0) { /* WHERE .LT. 0 MEANS ABSCISSAE COALESCED DURING THE QUADRATURE * STEP. ACUM IS GARBAGE DURING THE QUADRATURE STEP SO USE PACUM. */ sintc.acum = sintc.pacum; sintc.aacum = fabsf( sintc.acum ); sintc.epsmin = sintc.pepsmn; } else { sintnc.tps = sintc.acum - sintc.pacum; sintc.acum = sintc.pacum; sintnc.tps = (sintc.re*sintnc.tps)/(1.0e0 - sintc.re); if (fabsf( sintc.re - sintc.rep ) >= 0.125e0*sintc.re) { /* DO NOT TRUST THE ERROR ESTIMATE. */ sintc.err += fabsf( sintnc.tps ); } else { /* TRUST THE ERROR ESTIMATE, APPLY AITKEN ACCELERATION. */ sintc.acum -= sintnc.tps; } } if (sintnc.iprint > 1) sinto( 10, work ); goto L_2250; /* INTEGRATION WAS SUCCESSFUL ON LAST INTERVAL. CHECK IF THE * COMPLETE INTERVAL HAS BEEN INTEGRATED. * */ L_2235: if ((sintc.err*sintc.rndc) > sintc.epsmin) sintc.delta *= 1000.0e0; goto L_2240; L_2237: sintc.err = sintc.errc; L_2240: sintc.dischk = max( sintc.dischk - 1, 0 ); L_2250: sintc.search = 1; sintc.i = sintc.part; if (sintc.part != 1) { if (Step[1]*Step[2] > 0.0e0) sintc.i = 1; } Errt[sintc.i] += fmaxf( sintc.err, sintc.epsmin ); /* At one time, we had max(err,epsmin,edmeps*result(i)), where edmeps * is 0.5*r1mach(4), to account for round-off introduced during * accumulation. This effect should be negligible. */ Xcdobt[sintc.i] = Xcdobt[sintc.i] || sintc.err > powf(sintc.epsmin/ sintc.aacum,sintnc.relobt)*sintc.aacum; sintnc.tps = fmaxf( sintc.epso, fabsf( sintc.releps ) ); sintc.epsr = Errt[1] + Errt[2] + sintc.edue2a + sintc.edue2b; sintc.epsr = fmaxf( sintnc.tps - sintc.epsr, fmaxf( 0.1e0*sintc.epsr, 0.0e0 ) ); sintc.fail = FALSE; sintc.kmin = 0; Result[sintc.i] += sintc.acum; if (sintnc.iprint > 1) sinto( 11, work ); if (sintc.discf > 0) { if (sintc.part == 1) { sintc.discf = 0; } else { if (!sintc.iend) goto L_2350; } if (sintnc.iprint > 0) sinto( 13, work ); } if (!sintc.iend) goto L_2350; if (sintc.part != 1) goto L_2612; L_2320: *answer = Result[1]; if (sintnc.binit < sintnc.ainit) *answer = -*answer; Work[1] = Errt[1] + sintc.edue2a + sintc.edue2b; /* DO BLOCK */ if (Work[1] > sintc.epso) { Iflag[1] = -2; if (Xcdobt[1]) Iflag[1] = -3; if (sintnc.reltol == 0) goto L_2340; sintnc.tps = fabsf( sintc.releps ); sintc.releps = fabsf( *answer ) - Work[1]; sintc.releps = fmaxf( fabsf( Work[sintnc.reltol + 2] )*sintc.releps, sintec.enzer ); if (Work[1] > sintc.releps) { if (Work[1] > sintnc.tps) goto L_2340; if (Iflag[1] >= -2) goto L_2340; goto L_10; } } Iflag[1] = -1; /* END BLOCK */ L_2340: sintc.where = 0; L_2345: if (sintnc.nfindx != 0 && sintnc.kdim == 1) Iflag[sintnc.nfindx] = sintc.nfeval; return; /* THE INTEGRATION IS COMPLETE THROUGH THE PANEL ENDING AT BLOCAL. * */ L_2350: *fata = *fatb; *fatas = *fatbs; sintc.tlen += -sintc.absdif - sintc.absdif; if (labs( sintc.discf ) >= 2) { sintc.discf = isign( 1, sintc.discf ); *alocal = *blocal; sintc.eps = sintc.epsr*(fabsf( End[sintc.part] - *alocal )/ sintc.tlen); switch (IARITHIF(sintc.part - 1)) { case -1: goto L_40; case 0: goto L_40; case 1: goto L_50; } } sintc.kmax = maxk; sintnc.kk = sintc.k; if (sintc.did1) { if (sintc.k < sintc.kaimt) { if (sintc.err > sintec.ersqe6*sintc.epsmin) sintnc.kk += 1; } } else if (sintc.k < 6) { sintnc.kk = min( 6, sintc.k + sintc.nsub/2 ); } if (sintc.nsub != 0) { /* DO BLOCK */ if (fabsf( *blocal - Start[sintc.part] ) < fabsf( sintc.tend - Start[sintc.part] ) - sintc.delmin) { if (!sintc.havdif) goto L_2390; if (sintc.phisum < (float)( 10*(sintc.k - sintc.kaimt) )* sintc.phtsum) goto L_2380; if (sintc.phisum >= 16.0e0*sintc.phtsum) goto L_2390; if (sintc.phisum >= sintc.phtsum*powif(0.5e0,sintc.kaimt - sintc.k + 1)) goto L_2390; goto L_2380; } sintc.tend = End[sintc.part]; sintc.kmax = sintnc.kmaxf; /* END BLOCK */ L_2380: ; sintc.havdif = FALSE; sintns( NSRB ); } L_2390: ; *alocal = *blocal; sintc.errc = fmaxf( sintc.err, sintc.erri ); sintc.epnoiz = fminf( sintc.errc, 4.0e0*sintc.epnoiz ); if (!sintc.roundf) { if (sintc.rep > 0.3e0) sintc.kmin = sintc.k - 1; sintnc.s = 10.0e0*sintc.erri; if (sintnc.s <= sintc.epsmin && sintc.rndc > sintnc.fer) sintc.rndc = fmaxf( sintnc.fer, sintnc.s*sintc.rndc/sintc.epsmin ); } sintnc.s = fmaxf( sintc.rep, 1.0e-5 ); tp1 = 1.0e0; /* DO BLOCK */ L_2420: tp1 *= 1.125e0; sintnc.s *= 6.0e0; if (sintnc.s < 0.4e0) goto L_2420; /* END BLOCK * DO BLOCK */ if (((sintc.k >= max( sintc.kaimt - 1, 3 ) && sintc.k <= aimfor) && sintc.erri <= sintc.epsmin) && sintc.rep <= 0.1e0) { sintc.i = aimfor - sintc.k + 4; sintnc.kk = sintc.k - 1; } else { sintc.i = aimfor - sintnc.kk + 3; /* 1 .LE. I .LE. 7 */ if (sintc.i >= 3) { if ((sintc.i == 3 && sintc.erri >= sintc.epsmin) && sintc.rep > 0.3e0) goto L_2450; if (sintc.err == 0.0e0) goto L_2450; if (sintc.err > sintc.epsmin) goto L_2450; if (sintc.erri > 10.0e0*sintc.epsmin) goto L_2450; /* ADJUST STEP SIZE WHEN APPROACHING NOISE OR ROUND-OFF * LIMIT. */ sintc.delta *= tp1; goto L_2450; } } sintc.errc = fmaxf( sintc.errc, fabsf( sintc.ep ) ); sintc.err = sintc.errc; sintc.re = sintc.rep; if (sintc.rep > 1.0e0) sintc.delta *= 0.875e0; /* END BLOCK */ L_2450: ; sintnc.tps = sintc.eps; if (sintc.k >= sintnc.kk) sintnc.tps = fminf( sintec.esqeps*sintc.epsmin/sintc.rndc, sintc.eps ); alpha = 0.0e0; if (sintnc.tps != 0.0e0) alpha = fmaxf( sintc.errc*fminf( 1.0e0, sqrtf( sintc.re ) ), sintc.epsmin )/sintnc.tps; L_2480: sintc.delta *= Xstep[sintc.i]; sintnc.tp = fabsf( *alocal - Start[sintc.part] ); sintc.kaimt = aimfor; /* DO BLOCK */ if (sintc.search != 1 || sintnc.tp < fabsf( Worry[sintc.part] - Start[sintc.part] ) - sintc.delmin) { /* CHOOSE THE NEXT STEPSIZE BASED UPON THE ORDER OF THE CURRENT * FORMULA (SOME FUNCTION OF K), AND THE RATIO OF THE ERROR * COMMITTED TO THE ERROR REQUESTED ON THE CURRENT PANEL. * THE ERROR BEHAVES AS H**ORDER(K-1) WHERE H IS THE STEPSIZE. * THEN SET GAMMA = 1.125**ORDER(K-1), AND CHANGE THE ERROR * RATIO BY GAMMA AND THE STEPSIZE BY 1.125 UNTIL THE ERROR * RATIO IS GREATER THAN 2.0. * */ if (alpha != 1.0e0) { if (alpha < 1.0e0) { if (alpha == 0.0e0) goto L_2500; beta = 1.125e0; alpha *= Gamma[sintnc.kk - 1]; } else { alpha = 1.0e0/alpha; beta = 1.0e0/1.125e0; sintc.delta *= beta; } L_2490: ; alpha *= Gamma[sintnc.kk - 1]; if (alpha > 2.0e0) goto L_2500; sintc.delta *= beta; goto L_2490; } L_2500: ; if (sintc.search != 1) { sintc.dischk = max( sintc.dischk, 0 ); if (sintnc.iprint > 3) sinto( 12, work ); if (sintc.search == 6) goto L_2560; sintc.kmax = min( maxk - 1, sintc.kmax ); if (sintc.delta <= Z[2]*sintc.absdif || (sintc.k < 6 && sintc.delta < Z[3]*sintc.absdif)) goto L_2610; if (sintc.delta > Z[4]*sintc.absdif) sintc.kmax = min( sintc.kmax, 5 ); sintc.kaimt = min( sintc.kmax, sintc.kaimt ); if (sintc.k >= sintc.kmax) goto L_2560; sintc.delta = sintc.absdif + sintc.absdif; sintc.err = sintc.errc + sintc.extra; if (sintc.k >= 4) sintc.dischk = -1; if (sintc.err > sintc.eps) goto L_100; if (sintc.re > 1.0e0) goto L_100; switch (IARITHIF(konvrg - 6)) { case -1: goto L_100; case 0: goto L_100; case 1: goto L_2240; } } if (sintnc.kk < sintc.k) sintc.delta = fmaxf( sintc.delta, 2.0e0*tp1*Xstep[aimfor - sintc.k + 3]* sintc.absdif ); if (sintc.delta > sintc.absdif + sintc.absdif) sintc.kmin = sintc.k; if (sintnc.tp + sintc.delta < fabsf( End[sintc.part] - Start[sintc.part] )) { if (fminf( fabsf( sintc.tb ), sintnc.tp + sintc.delta ) >= fabsf( Worry[sintc.part] - Start[sintc.part] )) { sintc.dischk = -1; sintc.kmax = sintnc.kmaxf; if (sintnc.tp + 1.25e0*sintc.delta >= fabsf( sintc.tb )) goto L_2540; } } /* DO BLOCK */ if (sintnc.tp < fabsf( sintc.tend - Start[sintc.part] ) - sintc.delmin) { if (sintc.k >= sintc.kaimt) goto L_2560; if (!sintc.havdif) goto L_2560; if (sintc.phisum <= 10.0e0*sintc.phtsum) goto L_2560; if (sintc.nsub >= nsubmx) goto L_2560; *blocal = *alocal + signf( sintc.delta, sintc.tb ); sintns( NSIB ); goto L_2620; } sintc.tend = *alocal + signf( sintc.delta, sintc.tb ); sintc.kmax = sintnc.kmaxf; /* END BLOCK */ if ((sintc.did1 || sintc.part == 1) || sintc.fail) goto L_2610; /* AFTER THE FIRST SUCCESS IN PART 2, TREAT THE NEXT INTERVAL * LIKE A RESTART. ALSO, INCREASE THE STARTING STEPSIZE FOR * RESTARTING PART 1 IF POSSIBLE. */ sintc.kaimt = 3; sintc.kmax = sintnc.kmaxf; sintc.x1 = *alocal + signf( sintc.delta, sintc.tb ); sintc.x2 = *alocal; if (sintc.nsub != 0) { sintc.x1 = TDECR( sintc.x1 ); sintc.x2 = TDECR( sintc.x2 ); if (sintc.nsub != 2) { sintc.x1 = TDECR( sintc.x1 ); sintc.x2 = TDECR( sintc.x2 ); } } Step[1] = signf( fmaxf( fabsf( Step[1] ), fabsf( sintc.x1 - sintc.x2 ) ), Step[1] ); goto L_2610; } L_2540: ; sintc.dischk = -1; sintc.kmax = sintnc.kmaxf; sintc.delta = fmaxf( sintc.delta, fabsf( End[sintc.part] - *alocal ) ); if (sintc.delta > 3.0e0*sintc.absdif) sintc.dischk = -2; Worry[sintc.part] = End[sintc.part]; goto L_2610; /* END BLOCK */ L_2560: ; sintc.delta = fminf( sintc.absdif, sintc.delta ); Worry[sintc.part] = *blocal; sintc.kmax = sintnc.kmaxf; sintc.kaimt = min( max( 5, sintnc.kmaxf ), aimfor ); L_2610: *fatbs = FALSE; goto L_2620; L_2612: sintc.iend = FALSE; sintc.part = 1; sintc.nsub = 0; sintc.discf = 0; sintc.dischk = max( sintc.dischk, 1 ); Result[1] += Result[2]; Errt[1] += Errt[2]; Errt[2] = 0.0e0; Xcdobt[1] = Xcdobt[1] || Xcdobt[2]; *alocal = Start[1]; sintc.ta = sintc.tasave; sintc.tb = End[1] - sintc.ta; sintc.tlen = fabsf( End[1] - Start[1] ); sintc.delmin = sintsm( *alocal ); L_2614: if (sintc.nsub < sintc.nsubsv) { sintns( NSIA ); goto L_2614; } if (sintc.did1) sintc.delta = 0.0e0; sintc.delta = fmaxf( fabsf( Step[1] ), fmaxf( 10.0e0*sintc.delmin, sintc.delta ) ); Step[1] = signf( 1.0e0, signf( 1.0e0, sintc.tb )*(sintnc.binit - sintnc.ainit) ); *fatbs = FALSE; *fatas = sintc.fsaved; *fata = sintc.fsave; sintc.fsaved = FALSE; sintc.did1 = FALSE; L_2615: sintc.kaimt = 3; sintc.kmax = sintnc.kmaxf; L_2620: sintc.x2 = fabsf( End[sintc.part] - *alocal ); /* X2 IS THE LENGTH OF THE REST OF THE PART. */ sintc.init = FALSE; sintc.eps = sintc.epsr*sintc.x2/sintc.tlen; sintc.delmin = sintsm( *alocal ); if (sintc.delta < sintc.x2) { sintc.iend = FALSE; if (sintc.delta < sintc.delmin) { sintc.delta = sintc.delmin; *fatbs = FALSE; } /* DO NOT LEAVE A PIECE TOO SMALL. */ if (sintc.x2 - sintc.delta <= sintc.delmin) goto L_60; *blocal = *alocal + signf( sintc.delta, sintc.tb ); goto L_75; } /* REDUCE KAIMT IF DELTA IS VERY MUCH GREATER THAN THE REMAINING * PART OF THE MAJOR SUBDIVISION. * DO FOREVER */ L_2640: ; sintc.x2 *= 3.0e0; if (sintc.x2 > sintc.delta) goto L_60; sintc.kaimt -= 1; if (sintc.kaimt < 3) goto L_60; goto L_2640; /* END FOREVER * * ***** INTERNAL SUBROUTINES ******************************* * * COMPUTE MINIMUM RATIO OF STEPSIZE CHANGES. IF THE STEPSIZE * CHANGES MORE RAPIDLY THAN 0.3, ADD A POINT IN THE LARGEST * INTERVAL. * */ L_2670: sintc.where = 3; sintnc.kk = max( sintnc.inc2 + 1, 3 ); sintc.j = max( sintnc.inc, 3 ); sintc.eps = sintnc.s; goto L_2710; L_2680: sintc.where = 2; /* ABSCISSAE HAVE COALESCED. */ sintc.j = sintc.lendt; goto L_2700; L_2690: sintc.where = 1; sintc.j = sintc.lendt; if (jumps != 1) goto L_2700; if (sintc.search != 2 && sintc.search != 6) switch (IARITHIF(konvrg + konvrg - sintc.lendt)) { case -1: goto L_2700; case 0: goto L_2700; case 1: goto L_1190; } sintnc.kk = 4; sintc.j = 6; goto L_2710; L_2700: sintnc.kk = 3; L_2710: sintnc.tp = 0.3e0; tp1 = fabsf( Xt[sintnc.kk - 1] - Xt[sintnc.kk - 2] ); for (sintc.i = sintnc.kk; sintc.i <= sintc.j; sintc.i++) { sintnc.s = fabsf( Xt[sintc.i] - Xt[sintc.i - 1] ); if (sintnc.tp*sintnc.s >= tp1) { sintnc.tp = tp1/sintnc.s; sintc.l = sintc.i; } else if (sintnc.tp*tp1 > sintnc.s) { sintnc.tp = sintnc.s/tp1; sintc.l = sintc.i - 1; } tp1 = sintnc.s; } if (sintnc.tp == 0.3e0) switch (sintc.where) { case 1: goto L_910; case 2: goto L_1000; case 3: goto L_1160; } sintc.abscis = Xt[sintc.l - 1] + 0.5e0*(Xt[sintc.l] - Xt[sintc.l - 1]); if (sintc.where != 2) { if (sintc.abscis == Xt[sintc.l - 1]) goto L_1000; if (sintc.abscis == Xt[sintc.l]) goto L_1000; } if (jumps == 1) goto L_1520; if (sintc.j1 == 18) sintc.j1 = sintc.j2 + 1; if (sintc.j2 == 0) sintc.j2 = sintc.j1 - 1; goto L_1520; /* REQUEST A FUNCTION VALUE AT THE EDGE OF A PANEL. * */ L_2730: if (sintc.dischk <= 0) goto L_100; L_2740: switch (IARITHIF(sintc.endpts - 2)) { case -1: goto L_2750; case 0: goto L_2760; case 1: goto L_100; } L_2750: if (sintc.j1 == sintc.lendt) goto L_2760; sintc.where = 5; sintc.where2 = 1; sintc.endpts = 2; if (sintc.j2 == 1) sintc.endpts = 3; goto L_2810; L_2760: sintc.where = 5; sintc.where2 = 2; sintc.endpts = 3; if (sintc.j2 == 1) goto L_530; L_2810: sintc.abscis = Local[sintc.where2]; sintc.x = sintc.abscis; if (Fats[sintc.where2]) goto L_2830; goto L_2850; L_2820: Fat[sintc.where2] = sintnc.fncval; Fats[sintc.where2] = TRUE; Errat[sintc.where2] = sintnc.fer*fabsf( sintnc.fncval ); if (sintnc.fea != 0) Errat[sintc.where2] += fabsf( sintc.t*Work[sintnc.fea] ); if (Local[sintc.where2] == sintnc.ainit) { sintc.edue2a = sintnc.fncval*sintc.errina; sintc.eps = fmaxf( sintc.epsr - (sintc.edue2a + sintc.edue2b), sintc.edue2a + sintc.edue2b )*sintc.delta/sintc.tlen; } else if (Local[sintc.where2] == sintnc.binit) { sintc.edue2b = sintnc.fncval*sintc.errinb; sintc.eps = fmaxf( sintc.epsr - (sintc.edue2a + sintc.edue2b), sintc.edue2a + sintc.edue2b )*sintc.delta/sintc.tlen; } sintc.eps = fmaxf( sintc.eps, powf(sintc.epsmin/sintc.aacum,sintnc.relobt)* sintc.aacum ); L_2830: sintc.j = sintc.where2 + 2*(sintc.where - 4); switch (sintc.j) { case 1: goto L_280; case 2: goto L_290; case 3: goto L_350; case 4: goto L_350; } /* COMPUTE F USING FORWARD OR REVERSE COMMUNICATION. * */ L_2840: if (fabsf( sintc.abscis - Local[4 - sintc.where] ) > sintec.esmall) goto L_2850; sintc.where = -sintc.where; goto L_2180; L_2850: if (sintnc.nfmax <= 0) goto L_2860; if (sintc.nfeval < sintnc.nfmax) goto L_2860; sintc.where = 7; Iflag[1] = 5; goto L_2345; L_2860: sintc.nfeval += 1; sintc.t = 1.0e0; if (sintc.nsub != 0) { /* Put ABSCIS in user coordinates */ sintc.x1 = sintc.ta/sintc.tb; sintc.t = (sintc.abscis - sintc.ta)/sintc.tb; /*? ABSCIS=TA*(1.0e0+X1)+ABSCIS*(T-X1) */ sintc.abscis = sintc.ta + sintc.abscis*sintc.t + sintc.x1* (sintc.ta - sintc.abscis); if (sintc.nsub != 2) { sintc.f2 = (sintc.abscis - sintc.ta)/sintc.tb; /*? ABSCIS=TA*(1.0e0+X1)+ABSCIS*(F2-X1) */ sintc.abscis = sintc.ta + sintc.abscis*sintc.f2 + sintc.x1* (sintc.ta - sintc.abscis); sintc.t *= sintc.f2; sintc.t += sintc.t; } sintc.t += sintc.t; } if (sintc.abscis == *alocal) { /*## F2=SINTSM(ABSCIS)/EDELM2 *## ABSCIS=ABSCIS+SIGN(F2,TB) */ sintc.abscis += signf( fmaxf( sintec.enzer, sintec.emeps*fabsf( sintc.abscis ) ), sintc.tb ); sintc.x = sintc.abscis; } else if (sintc.abscis == *blocal) { /*## F2=SINTSM(ABSCIS)/EDELM2 *## ABSCIS=ABSCIS-SIGN(F2,TB) */ sintc.abscis -= signf( fmaxf( sintec.enzer, sintec.emeps*fabsf( sintc.abscis ) ), sintc.tb ); sintc.x = sintc.abscis; } Work[sintnc.kdim] = sintc.abscis; if (sintnc.revers != 0) return; sintf( answer, work, &Iflag[1] ); L_2900: sintnc.fncval = sintc.t**answer; switch (sintc.where) { case 1: goto L_90; case 2: goto L_150; case 3: goto L_160; case 4: goto L_2820; case 5: goto L_2820; case 6: goto L_350; case 7: goto L_2860; } #undef TINCR #undef TDECR } /* end of function */