/*Translated by FOR_C, v3.4.2 (-), on 07/09/115 at 08:31:20 */ /*FOR_C Options SET: ftn=u io=c no=p op=aimnv s=dbov str=l x=f - prototypes */ #include #include "fcrt.h" #include "sinto.h" #include #include /* PARAMETER translations */ #define IDFLT1 1024 #define IDINT1 (IDINT2*32) #define IDINT2 (IDFLT1*256) #define LDELMN 9 #define LDELTA 8 #define LDISCX 11 #define LEPS 167 #define LERRI 164 #define LINC 2 #define LISTOP 4 #define LK 2 #define LKAIMT 3 #define LKDIM 1 #define LLOC1 2 #define LLOC3 (LLOC1 + 2) #define LNSUB 4 #define LPART 5 #define LPHISU 48 #define LRNDC 161 #define LSEARC 6 #define LTPS 9 #define LTXTAA 1 #define LTXTAB 30 #define LTXTAC 62 #define LTXTAD 85 #define LTXTAE 106 #define LTXTAG 123 #define LTXTAH 157 #define LTXTAI 186 #define LTXTAJ 220 #define LTXTAK 239 #define LTXTAL 241 #define LTXTAM 264 #define LTXTAN 285 #define LTXTAO 304 #define LTXTAP 352 #define LTXTAQ 409 #define LTXTAR 439 #define LWHERE 7 #define LWORRY 66 #define LX 70 #define LXJ 10 #define LXT 73 #define MECONT 50 #define MEEMES 52 #define MEFDAT 25 #define MEFMAT 62 #define MEFVEC 61 #define MENTXT 23 #define MERET 51 #define MESA1 (IDINT1*LNSUB + (12 + LKDIM)*IDINT2 + IDFLT1*LLOC1 + LTXTAA) #define MESA10 (64*IDINT1 + IDFLT1*LRNDC + LTXTAN) #define MESA11 (64*IDINT1 + IDFLT1*LLOC3 + LTXTAO) #define MESA12 (IDINT1*(64 + LWHERE) + IDINT2*LNSUB + IDFLT1*LDELMN + LTXTAP) #define MESA13 (31*IDINT1) #define MESA14 (32*IDINT1 + IDFLT1*LDELTA + LTXTAQ) #define MESA15 (64*IDINT1 + IDFLT1*LDISCX + LTXTAR) #define MESA2 (IDFLT1*LDELTA + LTXTAB) #define MESA3 (IDINT1*LPART + IDINT2*LKAIMT + IDFLT1*LEPS + LTXTAC) #define MESA4 (29*IDINT1 + IDFLT1*LWORRY + LTXTAD) #define MESA5 (30*IDINT1 + LTXTAE) #define MESA6 ((32 + LK)*IDINT1 + (12 + LKDIM)*IDINT2 + IDFLT1*LERRI + LTXTAG) #define MESA7 (32*IDINT1 + LTXTAH) #define MESA8 (32*IDINT1 + IDFLT1*LX + LTXTAM) #define MESA9 (IDINT1*(96 + LNSUB) + IDFLT1*LDISCX + LTXTAA) #define MESL13 (-(15*IDINT1 + 3*IDINT2 + 25)) #define MESL3 (-(16*IDINT1 + 5*IDINT2 + 57)) #define MESL5 (-(IDINT2*LISTOP + IDFLT1*LXJ + LTXTAI)) #define MESL6 (-(8*IDINT1 + IDINT2*LINC + IDFLT1*LTPS + LTXTAL)) #define MESL8 (-(10*IDINT1 + 1*IDINT2 + 24)) #define MESL9 (-(11*IDINT1 + 6*IDINT2 + 68)) #define METEXT 53 /* 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*/ sinto( long jump, float work[]) { long int idat[5], if1, ii1, ii2, lidat, lmact, lmess, messf; float fdat[4]; static char mtxtaa[3][171]={"NSUB=$I on ($F, $F) KDIM=$I$B DEL=$(E12.5) DELMIN=$(E09.3) $BEPS=$G PART=$I AIM=$I$B DISCHK=$I WORRY=$F$B KDIM=$I$BTA=$F$E$(I2) $(E11.4) $G $G $G $G $G $J$E**** Reverse D", "irection ****$EISTOP=$I $I $I $I, XJ=$F XJP=$F$BXT$HFT$HPHI$HPHIT$E$#INC=$I INC2=$I E=$F $BX=$F F1=$F COUNT=$F$ERound-Off = $F. $BApparent non-integrable singularity n","ear $F. $BAbsiccae have coalesced. $ WHERE=$I, DELMIN=$F, NSUB=$I$EDELTA chosen by search = $F.$EDiscontinuity in ($F, $F).$EUsed $I function values -- the maximium.$E "}; static char mtxtab[1][81]={" K ERRI ERR EPSMIN EPS $ RE REPROD KDIM$E"}; static char mtxtac[1][33]={" PHISUM=$F PHTSUM=$F SEARCH=$I$E"}; static char mtxtad[1][68]={"Accept result $I = $F, ACCUM=$F, ERR=$(E10.4), EPSMIN=$G, KDIM=$I$E"}; static char mtxtae[1][38]={"Abscissae for dimensions $I to $I: $B"}; static char mtxtaf[1][59]={"Really on ($F, $F) KDIM=$I DEL=$(E12.5) TA=$(E10.3) $E"}; static char mtxtag[1][7]={"SINT$B"}; static long messl[13]={1,6,MESL3,7,MESL5,MESL6,9,MESL8,MESL9,12, 13,14,MESL13}; static long messa[15]={MESA1,MESA2,MESA3,MESA4,MESA5,MESA6,MESA7, MESA8,MESA9,MESA10,MESA11,MESA12,MESA13,MESA14,MESA15}; static long mact[26]={MENTXT,0,MEFDAT,0,METEXT,MENTXT,0,MEFDAT, 0,METEXT,MENTXT,0,MEFDAT,0,METEXT,MENTXT,0,MEFDAT,0,METEXT,MENTXT, 0,MEFDAT,0,METEXT,MERET}; static long mactma[7]={MEFMAT,17,0,4,LTXTAK,LTXTAJ,MERET}; static long mactar[6]={METEXT,MEFDAT,0,MEFVEC,0,MERET}; static long macth[2]={METEXT,MERET}; static long macter[5]={MEEMES,0,0,-1,MECONT}; float *const fnsav = (float*)&sintnc.ainit; float *const fsav = (float*)&sintc.aacum; long int *const insav = (long*)&sintnc.kdim; long int *const isav = (long*)&sintc.i; /* OFFSET Vectors w/subscript range: 1 to dimension */ float *const Fdat = &fdat[0] - 1; float *const Fnsav = &fnsav[0] - 1; float *const Fsav = &fsav[0] - 1; long *const Idat = &idat[0] - 1; long *const Insav = &insav[0] - 1; long *const Isav = &isav[0] - 1; float *const Local = &sintc.local[0] - 1; long *const Mact = &mact[0] - 1; long *const Mactar = &mactar[0] - 1; long *const Macter = &macter[0] - 1; long *const Macth = &macth[0] - 1; long *const Mactma = &mactma[0] - 1; long *const Messa = &messa[0] - 1; long *const Messl = &messl[0] - 1; double *const Result = &sintc.result[0] - 1; float *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. *>> 2009-07-15 SINTO Krogh Fixed incorrect outputs of KDIM *>> 2008-01-17 SINTO Krogh Updated generated error message text. *>> 2000-12-01 SINTO Krogh Removed unused parameters LABSC & MASA16. *>> 1996-03-31 SINTO Krogh Removed unused variable in common. *>> 1995-11-20 SINTO Krogh Converted from SFTRAN to Fortran 77. *>> 1994-11-17 SINTO Krogh Rearranged parameter statments. *>> 1994-11-14 SINTO Krogh Declared all vars. *>> 1994-10-19 SINTO Krogh Changes to use M77CON *>> 1994-07-07 SINTO Snyder set up for CHGTYP. *>> 1994-06-13 SINTO Krogh -- Fixed value of LLOC3 "+2" not "+3". *>> 1993-05-18 SINTO Krogh -- Changed "END" to "END PROGRAM" *>> 1993-04-29 SINTO Krogh Additions for Conversion to C. *>> 1993-04-13 SINTO Krogh Minor changes for new MESS. *>> 1992-05-28 SINTO Krogh Corrected minor problem in error message. *>> 1992-04-08 SINTO Krogh unused label 210, FDAT in MESS calls removed *>> 1992-03-03 SINTO Krogh converted to use message processor. *>> 1991-09-20 SINTO Krogh converted '(1)' dimensioning to '(*)'. *>> 1988-11-17 SINTO Snyder Remove unnecessary specializer comments. *>> 1988-06-07 SINTO Snyder Correct format statement 20. *>> 1987-11-19 SINTO Snyder Initial code. * * PRINT SOMETHING FOR SINTA. * *--S replaces "?": ?INT, ?INTA, ?intc, ?intec, ?INTNC, ?INTO, ?MESS * * ***** DESCRIPTION OF VARIABLES WITH SPECIAL USE HERE ********* * * AACUM First floating point variable in saved common area. * ACUM Double precision result accumulated. * AINIT First floating point variable in the unsaved common area. * DISCHK Nonzero if checking for a discontinuity. This and WORRY are * only printed if DISCHK is not zero. * SMESS Used when message contains floating point. * FDAT Temporary storage for floating point variables to be printed. * FNSAV Array equivalenced to AINIT, access to unsaved floating point. * FSAV Array equivalenced to AACUM, access to saved floating point. * I First integer in the saved common area. * IDAT Temporary storage for integer variables to be printed. * IDFLT1 Parameter, see description of MESSA. * IDINT1 Parameter, see description of MESSA. * IDINT2 Parameter, see description of MESSA. * IF1 Pointer to floating point variable, see MESSA. * II1 First integer extracted from MESSF, see MESSA. * II2 Second integer extracted from MESSF, see MESSA. * INSAV Array equivalenced to KDIM for access to unsaved integers. * ISAV Array equivalenced to I for access to saved integers. * JUMP Input variable defining what is to be printed. * = 1 Panel boundaries * = 2 Header (if needed) and K, ERRI, ERR, EPSMIN, EPS, RE, ... * = 3 No action. * = 4 Note that direction of accumulation has been reversed. * = 5 x's, f's, difference lines, etc. * = 6 Estimated errors during a search. * = 7 Panel boundaries after a disconinuity found. * = 8 New round off level after noise detected. * = 9 Indicate a nonintegrable singularity. * =10 Note that abscissae have coalesced. * =11 Data for an accepted answer. * =12 Step size from an initial interval search. * =13 Message that there appears to be a discontinuity. * KDIM First integer variable in the unsaved common area. * L?? Many names starting with L identify the locations of variables * the common blocks. The letters following the L serve to * identify which variable in the common blocks. All of these * variables are defined in a block of parameter statements * separated from others by header and ending comments. * LENDT Number of x's, f,s and differences to be printed. * LIDAT Current location for saving integers in IDAT. * LMACT Current base location for saving data in MACT. * LMESS Value taken from MESSL and updated depending on MESSF, see * description of MESSA below. If > 0, provides an index * into MESSA for the action. After getting actions, if this * is > 1, it gives the next value of LMESS; else data is printed * and a return is made with printout of exterior abscissae * following the other data if LMESS is 1. If LMESS < 0, -LMESS * is packed data like that in MESSA, except that in this case * II1 gives the next value to be assigned to LMESS, and both the * integer and floating point data is to be printed from the * unsaved common block. * LTXTA? All names of this type were generated by running the data in * SINTO.ERR through PMESS. These are all parameters that define * where various text starts in MTXTAx. * MACT Array used to store actions for calls to MESS. * MACTAR As for MACT except for printing a vector. * MACTH As for MACT except only need pointer to one text entry. * MACTMA As for MACT except for printing the array containing x's, f's, * and differences. * MECONT Parameter defined in MESS, means return, print is to continue. * MEFDAT Parameter defined in MESS, means set index for next floating * point item to print. * MEFMAT Parameter defined in MESS, means print a matrix. * MEFVEC Parameter defined in MESS, means print a vector. * MENTXT Parameter defined in MESS, as for MEFDAT execpt for text. * MERET Parameter defined in MESS, means print buffers and return. * MESAdd Where dd is one or two digits. These define parameters that * are used in the data statement for MESSA, see comments there. * MESLd As for the above, execpt used in assigning value to MESSL, when */ /* those value are < 0. * MESS Message routine, when no floating point is to be printed. * MESSA An array containing packed integers that define the actions to * be taken. Entries in MESSA have the form IDINT1*II1 + * IDINT2*II2 + IDFLT1*IF1 + t, where in every case the first * multiplier in a product is a parameter which is a power of 2, * and the second defines an action as follows: * II1 first integer to be printed. In most cases, LMESS increases * by 1, but if this is > 31, LMESS is set to II1 / 32 - 1 and * II1 is replaced by mod(II1,32). If the result is: * 31 Print out of final results (some vars. are always double * precision). * 30 Check on NDIM to see if KDIM should print. * 29 Check on DISCHK to decide if it and WORRY are to print. * >12 Subtract 12, and use to get an integer from unsaved common. * <12 Use to get an integer from the saved common area. * II2 Second integer. Treated like the first, except always < 25. * IF1 Used to get a floating point number from the saved common * block. * t Location of text for message in MTXTAA. * MESSF Obtained from MESSA or -MESSL(LMESS), see above. * MESSL Maps value of JUMP into actions. If > 0, gives an index into * MESSA, else the negative defines the value of MESSF. Also see * LMESS above. * METEXT Parameter defined in MESS, means print as defined by MTXTAx. * MTXTAx Character arrays containing text and instructions for the * printing of messages by MESS. * NDIM Number of dimensions in the total integral. * NEEDH .TRUE. if heading needed when JUMP is 2, else is .FALSE. * SMESS Used when message contains single precision floating point. * WORK Array passed in containing exterior abscissae. * * ***** PROGRAM VARIABLES ********************************** * */ /* ***** COMMON VARIABLES *********************************** * * 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 ***************************** * */ /* ***** Statements for Processing Messages ********************** * */ /* Parameters defining locations in the common blocks. */ /* End of parameters defining locations in the common blocks. * * ********* Error message text *************** *[Last 2 letters of Param. name] [Text generating message.] *AA NSUB=$I on ($F, $F) KDIM=$I$B *AB DEL=$(E12.5) DELMIN=$(E09.3) $B *AC EPS=$G PART=$I AIM=$I$B *AD DISCHK=$I WORRY=$F$B *AE KDIM=$I$B *AF TA=$F$E *AG $(I2) $(E11.4) $G $G $G $G $G $J$E *AH **** Reverse Direction ****$E *AI ISTOP=$I $I $I $I, XJ=$F XJP=$F$B *AJ XT$HFT$HPHI$HPHIT$E *AK $# *AL INC=$I INC2=$I E=$F $B *AM X=$F F1=$F COUNT=$F$E *AN Round-Off = $F. $B *AO Apparent non-integrable singularity near $F. $B *AP Absiccae have coalesced. WHERE=$I, DELMIN=$F, NSUB=$I$E *AQ DELTA chosen by search = $F.$E *AR Discontinuity in ($F, $F).$E *AS Used $I function values -- the maximium.$E * $ *AT K ERRI ERR EPSMIN EPS $C * RE REPROD KDIM$E * $ *AU PHISUM=$F PHTSUM=$F SEARCH=$I$E * $ *AV Accept result $I = $F, ACCUM=$F, ERR=$(E10.4), EPSMIN=$G,$C * KDIM=$I$E * $ *AW Abscissae for dimensions $I to $I: $B * $ *AX Really on ($F, $F) KDIM=$I DEL=$(E12.5) TA=$(E10.3) $E * $ *AY SINT$B */ /* **** End of automatically generated text */ /* JUMP = 1 2 3 4 5 6 7 8 9 */ /* JUMP = 10 11 12 13 */ /* ***** EXECUTABLE STATEMENTS ****************************** * */ if (jump == 2) { if (sintnc.needh) mess( macth, (char*)mtxtab,81, idat ); sintnc.needh = FALSE; } else { sintnc.needh = TRUE; } lmact = 2; lidat = 1; lmess = Messl[jump]; if (lmess < 0) { /* Output data from the unsaved common block */ messf = -lmess; goto L_20; } L_10: messf = Messa[lmess]; L_20: ii1 = messf/IDINT1; messf += -ii1*IDINT1; ii2 = messf/IDINT2; messf += -ii2*IDINT2; if1 = messf/IDFLT1; Mact[lmact] = messf - if1*IDFLT1; Mact[lmact + 2] = if1; if (lmess <= 0) { lmess = ii1; if (if1 == 0) { if (lmess == 0) return; Macter[2] = Mact[lmact]; Macter[3] = ii2; mess( macter, (char*)mtxtag,7, idat ); sintnc.ndim += 100; goto L_10; } Mact[6] = MECONT; smess( mact, (char*)mtxtaa,171, &Insav[ii2], fnsav ); Mact[6] = MENTXT; if (if1 != LXJ) goto L_10; smess( macth, (char*)mtxtac,33, &Isav[LSEARC], &Fsav[LPHISU] ); Mactma[3] = sintc.lendt; smess( mactma, (char*)mtxtaa,171, idat, &Fsav[LXT] ); return; } else { if (ii1 >= 32) { lmess = ii1/32; ii1 += -32*lmess; lmess -= 1; } else { lmess += 1; } if (ii1 != 0) { if (ii1 > 12) { if (ii1 >= 28) { switch (ii1 - 28) { case 1: goto L_100; case 2: goto L_110; case 3: goto L_120; } lmess = 1; goto L_220; L_100: if (sintc.dischk == 0) goto L_10; Idat[lidat] = sintc.dischk; Mact[lmact + 2] = if1 + sintc.part - 1; goto L_200; L_110: lmess = 1; if (sintc.nsub != 0) lmess = -1; if (sintnc.ndim == 1) goto L_220; Idat[lidat] = sintnc.kdim; goto L_200; /* Take care of stuff that is double precision in single precision code. */ L_120: Fdat[1] = Result[sintc.i]; Fdat[2] = sintc.acum; Fdat[3] = sintc.err; Fdat[4] = sintc.epsmin; Idat[1] = sintc.i; Idat[2] = sintnc.kdim; smess( macth, (char*)mtxtad,68, idat, fdat ); goto L_230; } else { Idat[lidat] = Insav[ii1 - 12]; } } else { Idat[lidat] = Isav[ii1]; } if (ii2 != 0) { lidat += 1; if (ii2 > 12) { Idat[lidat] = Insav[ii2 - 12]; } else { Idat[lidat] = Isav[ii2]; } } lidat += 1; } } L_200: lmact += 5; if (lmess > 1) goto L_10; L_220: Mact[lmact - 1] = MERET; if (sintnc.ndim > 100) Mact[lmact - 1] = MECONT; smess( mact, (char*)mtxtaa,171, idat, fsav ); Mact[lmact - 1] = MEFDAT; if (lmess == 0) return; if (lmess < 0) { sintc.abscis = fabsf( Local[4] - Local[3] ); Idat[1] = sintnc.kdim; smess( macth, (char*)mtxtaf,59, idat, &Fsav[LLOC3] ); } L_230: if (sintnc.ndim <= sintnc.kdim) return; lmact = 1; if (sintnc.ndim > 100) { sintnc.ndim -= 100; if (sintnc.ndim <= sintnc.kdim) lmact = 2; } Idat[2] = sintnc.ndim; Idat[1] = sintnc.kdim + 1; Mactar[3] = Idat[1]; Mactar[5] = -sintnc.ndim; smess( &Mactar[lmact], (char*)mtxtae,38, idat, work ); return; } /* end of function */