/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:32:03 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "dint1.h" #include /* COMMON translations */ struct t_dintnc { double 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; } dintnc; struct t_dintc { double acum, pacum, result[2], 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]; } dintc; struct t_dintec { double emeps, eepsm8, edelm2, edelm3, esqeps, ersqep, ersqe6, eminf, esmall, enzer, edelm1, eninf; } dintec; /* end of COMMON translations */ void /*FUNCTION*/ dint1( double a, double b, double *answer, double work[], long iopt[]) { /* OFFSET Vectors w/subscript range: 1 to dimension */ long *const Iopt = &iopt[0] - 1; double *const Work = &work[0] - 1; /* end of OFFSET VECTORS */ /* Copyright (c) 1996 California Institute of Technology, Pasadena, CA. * ALL RIGHTS RESERVED. * Based on Government Sponsored Research NAS7-03001. *>> 1996-03-31 DINT1 Krogh Removed unused variable in common. *>> 1995-11-20 DINT1 Krogh Converted from SFTRAN to Fortran 77. *>> 1994-10-19 DINT1 Krogh Changes to use M77CON *>> 1994-07-07 DINT1 Snyder set up for CHGTYP. *>> 1993-05-18 DINT1 Krogh -- Changed "END" to "END PROGRAM" *>> 1991-09-20 DINT1 Krogh converted '(1)' dimensioning to '(*)'. *>> 1987-11-19 DINT1 Snyder Initial code. * *--D replaces "?": ?INT1, ?INTA, ?intc, ?intec, ?INTF, ?INTOP, ?INTNC * ****************************************************************** * * THIS SUBROUTINE ATTEMPTS TO CALCULATE THE INTEGRAL OF A FUNCTION * OVER THE INTERVAL A TO B TO WITHIN A TOLERANCE LEVEL SPECIFIED * BY THE USER. THE FUNCTION IS PROVIDED BY THE USER VIA A * SUBPROGRAM REFERENCED BY CALL DINTF (F,X,IFLAG), OR BY REVERSE * COMMUNICATION. IN THE CASE OF ONE-DIMENSIONAL QUADRATURE, IFLAG * IS ALWAYS ZERO. ALL ABSCISSAE ARE WITHIN THE INTERVAL OF * INTEGRATION. * * 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 INTERVAL OF * INTEGRATION, AND NON-INTEGRABLE SINGULARITIES NEAR THIS INTERVAL * (INCLUDING COMPLEX VALUES). THE EFFICIENCY AND RELIABILITY * OF INTEGRATING SUCH FUNCTIONS MAY BE GREATLY IMPROVED BY MANUALLY * SUBDIVIDING THE INTERVAL OF INTEGRATION AT THE DISCONTINUITY, * SINGULARITY OR CLOSEST POINT TO A POLE, AND SUMMING THE ANSWERS. * A CHANGE OF VARIABLE TO ELIMINATE OR REDUCE THE STRENGTH OF THE * SINGULARITY WILL SIGNIFICANTLY IMPROVE PERFORMANCE. * * ***** FORMAL ARGUMENTS ************************************* * * A LOWER LIMIT OF INTERVAL OF INTEGRATION. */ /* B UPPER LIMIT OF INTERVAL OF INTEGRATION. */ /* ANSWER IS THE ESTIMATE OF THE INTEGRAL. WHEN REVERSE COMMUNICATION * IS SPECIFIED IT IS USED TO PASS FUNCTION VALUES FROM * THE USER PROGRAM TO DINTA. */ /* WORK IS A WORKING STORAGE VECTOR USED BY OPTIONS (SEE IOPT). * WHEN THE INTEGRATION IS COMPLETE, WORK(1) CONTAINS THE * ESTIMATED ABSOLUTE ERROR. WHEN REVERSE COMMUNICATION * IS SPECIFIED, WORK(1) IS USED TO PASS ABSCISSAE FROM * DINTA TO THE USER PROGRAM. */ /* IOPT IS A VECTOR OF INTEGERS USED TO RETURN THE STATUS AND TO * SPECIFY OPTIONS. */ /* IOPT(1) IS USED TO RETURN A STATUS INDICATOR TO THE USER. * VALUES OF THE STATUS INDICATOR ARE * -1 - NORMAL TERMINATION, EITHER THE ABSOLUTE OR THE RELATIVE * ERROR CRITERIA ARE SATISFIED. * -2 - NORMAL TERMINATION, NEITHER THE ABSOLUTE NOR THE * 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. * +5 - TOO MANY FUNCTION VALUES NEEDED. * +6 - APPARENT NON-INTEGRABLE SINGULARITY NEAR ANSWER. * ENTRIES IN IOPT STARTING AT IOPT(2) ARE DESCRIBED BELOW. * IOPT(I) ENTRY IN IOPT(I) MEANS * 0 NO MORE OPTIONS, BEGIN INTEGRATION. * 1 IOPT(I+1) CONTAINS THE UNIT NUMBER FOR OUTPUT. * 2 IOPT(I+1) IS DIAGNOSTIC PRINT LEVEL * 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). * 3 WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE TOLERANCE, * WORK(IOPT(I+1)+1) CONTAINS TOLERANCE RELATIVE TO * THE LOCALLY OBTAINABLE PRECISION, AND * WORK(IOPT(I+1)+2) CONTAINS TOLERANCE RELATIVE TO * THE VALUE OF THE INTEGRAL. THE TOLERANCE RELATIVE * TO THE LOCALLY OBTAINABLE PRECISION IS SPECIFIED AS * THE FRACTION OF LOCALLY OBTAINABLE DIGITS THAT ARE * PERMITTED TO BE WRONG, AND IS INTERNALLY BOUNDED * BETWEEN 0.0 AND 1.0. IF ANY OF THE ERROR CONTROL * CRITERIA ARE SATISFIED, THE ANSWER IS ACCEPTED. IF * THIS OPTION IS NOT USED, THE ABSOLUTE AND RELATIVE * TOLERANCES ARE SET TO ZERO, AND THE TOLERANCE * RELATIVE TO THE LOCALLY OBTAINABLE PRECISION IS SET * TO 0.25. * 4 WORK(IOPT(I+1)) CONTAINS ABSOLUTE ERROR COMMITTED * COMPUTING F. WORK(IOPT(I+1)) MAY BE CHANGED DURING * THE INTEGRATION. * 5 WORK(IOPT(I+1)) CONTAINS RELATIVE ERROR EXPECTED TO * BE COMMITTED COMPUTING F. CHANGES TO THIS VALUE * DURING INTEGRATION WILL NOT BE DETECTED. * 6 USE REVERSE COMMUNICATION - * CALL DINT1 (A,B,ANSWER,WORK,IOPT) * 1 CALL DINTA (ANSWER,WORK,IFLAG) * IF (IFLAG.NE.0) GO TO 2 * ANSWER=F(WORK(1)) * GO TO 1 * 2 CONTINUE * IFLAG AT THIS POINT CONTAINS THE SAME VALUE AS * THE STATUS FLAG WOULD (IF SELECTED) FOR THE * FORWARD COMMUNICATION METHOD. * 7 SET MINIMUM INDEX OF QUADRATURE FORMULA TO USE * BEFORE SUBDIVISION TO IOPT(I+1). * 8 NO EFFECT. MAY BE USED TO PASS INFORMATION INTO * DINTF. INCREMENT I BY 2. * 9 IOPT(I+1) IS THE MAXIMUM NUMBER OF FUNCTION VALUES * TO USE TO INTEGRATE THE ENTIRE PROBLEM. IF * IOPT(I+1) .LE. 0, THE NUMBER OF FUNCTION VALUES * IS NOT CONTROLLED. * 10 IOPT(I+1) IS USED TO RETURN THE NUMBER OF FUNCTION * VALUES USED TO INTEGRATE THE ENTIRE PROBLEM. * 11 WORK(IABS(IOPT(I+1))) IS THE LOCATION OF A * SINGULARITY OR DISCONTINUITY. IF THE LOCATION IS * INSIDE THE INTERVAL, THE INTERVAL WILL BE * IMMEDIATELY SUBDIVIDED. IF IOPT(I+1).GT.0, A * T**2 SUBSTITUTION BASED AT WORK(...) WILL BE USED. * IF IOPT(I+1) .LT. 0 A T**4 SUBSTITUTION WILL BE * USED. IF IOPT(I+1) .EQ. 0, NO SUBSTITUTION WILL BE * USED. * 12 WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE ERROR IN THE * LOWER LIMIT OF THE CURRENT DIMENSION. THE ERROR IN * THE UPPER LIMIT IS IN WORK(IOPT(I+1)+1). * 13 IS USED ONLY FOR MULTIDIMENSIONAL PROBLEMS. * * ALL OPTIONS ARE SET TO THEIR NOMINAL VALUES BEFORE THE OPTION * VECTOR IS PROCESSED. * * ***** EXTERNAL REFERENCES ********************************* * * DINTA TO PERFORM THE INTEGRATION. * DINTOP TO SET OPTIONS. * * ***** COMMON STORAGE ************************************** * * Common blocks are: DINTNC, DINTC, and DINTEC * COMMON /DINTNC/ CONTAINS VARIABLES NOT SEPARATELY SAVED FOR * EACH DIMENSION OF A MULTIPLE QUADRATURE. COMMON /DINTC/ * 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 /DINTC/. 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 /DINTNC/ VARIABLES. * */ /* DECLARATIONS OF COMMON /DINTC/ 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 DINTC. * 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. */ /* ***** PROCEDURES ****************************************** * */ dintc.where = 0; dintnc.ndim = 1; dintnc.kdim = 1; dintc.nfeval = 0; dintop( iopt, work ); if (Iopt[1] == 0) { /* CALL DINTA TO DO THE INTEGRATION. * */ dintnc.ainit = a; dintnc.binit = b; if (dintnc.revers == 0) dinta( answer, work, iopt ); } return; } /* end of function */