1      DOUBLE PRECISION FUNCTION DINTSM (SXMIN)
2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
3c ALL RIGHTS RESERVED.
4c Based on Government Sponsored Research NAS7-03001.
5C>> 1996-03-31 DINTSM  Krogh  Removed unused variable in common.
6c>> 1995-11-20 DINTSM  Krogh  Converted from SFTRAN to Fortran 77.
7c>> 1994-10-19 DINTSM  Krogh  Changes to use M77CON
8c>> 1994-07-07 DINTSM  Snyder set up for CHGTYP.
9c>> 1994-07-05 DINTSM  Snyder  Corrected calculation
10C>> 1993-05-18 DINTSM  Krogh -- Changed "END" to "END PROGRAM"
11C>> 1987-11-19 DINTSM  Snyder  Initial code.
12C
13c--D replaces "?": ?intc, ?intec, ?intnc, ?INTSM
14c
15C     CALCULATE THE MINIMUM STEPSIZE TO USE IF ALOCAL WERE SET EQUAL TO
16C     SXMIN.
17C
18C     WRITE X = TA + (T-TA)**2/TB.  IF WE LET X2 - X1 BE THE SMALLEST
19C     ALLOWED STEP AT X1, SAY SMIN, THEN
20C     TB*(X2-X1) = TB*SMIN = (T2-TA)**2 - (T1-TA)**2, OR
21C     TB*SMIN = (T2-T1)*(T2-T1+2*(T1-TA)).  SOLVING FOR T2-T1 PROVIDES
22C     THE EXPRESSIONS IN THE CODE BELOW.  THE ANALYSIS PROCEEDS
23C     SIMILARLY WHEN X = TA + (T-TA)**4/TB**3.
24C
25      DOUBLE PRECISION SXMIN
26C
27C     *****     LOCAL VARIABLES     ************************************
28C
29C SG      IS A TEMPORARY VARIABLE
30      DOUBLE PRECISION SG
31C SMIN    IS THE VALUE THAT WILL BE RETURNED AS THE MINIMUM STEPSIZE.
32      DOUBLE PRECISION SMIN
33C SOLVE   IS AN ARITHMETIC STATEMENT FUNCTION DEFINED BELOW.
34      DOUBLE PRECISION SOLVE
35C SQRTTB  IS SQRT(ABS(TB))
36      DOUBLE PRECISION SQRTTB
37C SX      IS A LOCAL COPY OF SXMIN.
38      DOUBLE PRECISION SX
39C TDECR   IS AN ARITHMETIC STATEMENT FUNCTION DEFINED BELOW.
40      DOUBLE PRECISION TDECR
41C
42C     *****     COMMON VARIABLES     ***********************************
43C
44C     COMMON /DINTNC/ CONTAINS VARIABLES NOT SEPARATELY SAVED FOR
45C     EACH DIMENSION OF A MULTIPLE QUADRATURE.  COMMON /DINTC/
46C     CONTAINS VARIABLES THAT MUST BE SAVED FOR EACH DIMENSION OF THE
47C     QUADRATURE.  THE VARIABLES IN EACH COMMON BLOCK ARE STORED IN THE
48C     ORDER - ALWAYS DOUBLE, DOUBLE IF DOUBLE PRECISION PROGRAM, DOUBLE
49C     IF DOUBLE PRECISION PROGRAM AND EXPONENT RANGE OF DOUBLE AND
50C     SINGLE VERY DIFFERENT, SINGLE, INTEGER, LOGICAL.  A PAD OF LOGICAL
51C     VARIABLES IS INCLUDED AT THE END OF /DINTC/.  THE DIMENSION OF
52C     THE PAD MAY NEED TO BE VARIED SO THAT NO VARIABLES BEYOND THE END
53C     OF THE COMMON BLOCK ARE ALTERED.
54C
55C     DECLARATIONS OF COMMON /DINTNC/ VARIABLES.
56C
57      DOUBLE PRECISION AINIT, BINIT, FNCVAL, S, TP
58      DOUBLE PRECISION FER, FER1, RELOBT, TPS, XJ, XJP
59      INTEGER     FEA,       FEA1,      INC,       INC2,      IPRINT,
60     1 ISTOP(2,2),JPRINT,    KDIM,      KK,        KMAXF,     NDIM,
61     2 NFINDX,    NFMAX,     NFMAXM,    RELTOL,    REVERM,    REVERS,
62     3 WHEREM
63      LOGICAL NEEDH
64C
65C     DECLARATIONS OF COMMON /DINTC/ VARIABLES.
66C
67c--D Next line special: S => D, X => Q, D => D, P => D
68      DOUBLE PRECISION ACUM, PACUM, RESULT(2)
69C     139 $.TYPE.$ VARIABLES
70      DOUBLE PRECISION
71     1 AACUM,     ABSCIS,    DELMIN,    DELTA,     DIFF,      DISCX(2),
72     2 END(2),    ERRINA,    ERRINB,    FAT(2),    FSAVE,
73     3 FUNCT(24), F1,        F2,        LOCAL(4),  PAACUM,    PF1,
74     4 PF2,       PHISUM,    PHTSUM,    PX,        SPACE(6),
75     5 STEP(2),   START(2),  SUM,       T,         TA,        TASAVE,
76     6 TB,        TEND,      WORRY(2),  X,         X1,
77     7 X2,        XT(17),    FT(17),    PHI(34)
78c Note XT, FT, and PHI above are last, because they must be in adjacent
79c locations in DINTC.
80C     30 $DSTYP$ VARIABLES
81      DOUBLE PRECISION
82     1 ABSDIF,    COUNT,     EDUE2A,    EDUE2B,    EP,        EPNOIZ,
83     2 EPS,       EPSMAX,    EPSMIN,    EPSO,      EPSR,      EPSS,
84     3 ERR,       ERRAT(2),  ERRC,      ERRF,      ERRI,      ERRT(2),
85     4 ESOLD,     EXTRA,     PEPSMN,    RE,        RELEPS,    REP,
86     5 REPROD,    RNDC,      TLEN,      XJUMP
87C     29 INTEGER VARIABLES
88      INTEGER     DISCF,     DISCHK,    ENDPTS,    I,         INEW,
89     1 IOLD,      IP,        IXKDIM,    J,         J1,        J1OLD,
90     2 J2,        J2OLD,     K,         KAIMT,     KMAX,      KMIN,
91     3 L,         LENDT,     NFEVAL,    NFJUMP,    NSUB,      NSUBSV,
92     4 NXKDIM,    PART,      SEARCH,    TALOC,     WHERE,     WHERE2
93C     11 TO 18 LOGICALS (7 ARE PADDING).
94      LOGICAL     DID1,      FAIL,      FATS(2),   FSAVED,    HAVDIF,
95     1 IEND,      INIT,      ROUNDF,    XCDOBT(2), PAD(7)
96C
97C     THE COMMON BLOCKS.
98C
99      COMMON /DINTNC/
100c        1       2       3     4        5       6       7        8
101     W AINIT,  BINIT,  FNCVAL, S,      TP,     FER,    FER1,   RELOBT,
102c       9      10       11      12      13       1       2        3
103     X TPS,    XJ,     XJP,    FEA,    FEA1,   KDIM,    INC,    INC2,
104c     4 (2,2)    8       9     10       11      12       13      14
105     Y ISTOP,  JPRINT, IPRINT, KK,     KMAXF,  NDIM,   NFINDX, NFMAX,
106c        15     16       17      18      19      20
107     Z NFMAXM, RELTOL, REVERM, REVERS, WHEREM, NEEDH
108      COMMON /DINTC/
109     1 ACUM,   PACUM,  RESULT
110      COMMON /DINTC/
111c        1     2 (4)     6      7        8       9      10     11 (2)
112     1 AACUM,  LOCAL,  ABSCIS, TA,     DELTA,  DELMIN, DIFF,   DISCX,
113c     13 (2)     15      16    17 (2)   19     20 (24) 44
114     2 END,    ERRINA, ERRINB, FAT,    FSAVE,  FUNCT,  F2,
115c       45      46     47       48      49     50      51 (6)
116     3 PAACUM, PF1,    PF2,    PHISUM, PHTSUM, PX,     SPACE,
117c      57 (2)  59 (2)   61     62        63    64       65
118     4 STEP,   START,  SUM,    T,      TASAVE, TB,     TEND,
119c      66 (2)  68      69      70      71       72
120     5 WORRY,  X1,     X2,     X,      F1,     COUNT,
121c      73 (17) 90 (17) 107 (34)
122     6 XT,     FT,     PHI
123      COMMON /DINTC/
124c       141     142    143     144      145     146
125     1 ABSDIF, EDUE2A, EDUE2B, EP,     EPNOIZ, EPSMAX,
126c       147     148     149    150 (2)  152     153
127     2 EPSO,   EPSR,   EPSS,   ERRAT,  ERRC,   ERRF,
128c     154 (2)   156     157     158     159    160
129     3 ERRT,   ESOLD,  EXTRA,  PEPSMN, RELEPS, REP,
130c       161     162     163
131     4 RNDC,   TLEN,   XJUMP,
132c       164    165      166    167    168       169
133     5 ERRI,   ERR,    EPSMIN, EPS,    RE,     REPROD
134      COMMON /DINTC/
135c       170     171     172
136     1 DISCF,  DISCHK, ENDPTS, INEW,   IOLD,   IP,     IXKDIM,
137     2 J,      J1,     J1OLD,  J2,     J2OLD,  KMAX,   KMIN,
138     3 L,      LENDT,  NFJUMP, NSUBSV, NXKDIM, TALOC,  WHERE2,
139c      1       2          3      4       5         6      7       8
140     4 I,      K,      KAIMT,  NSUB,   PART,   SEARCH, WHERE, NFEVAL
141      COMMON /DINTC/
142     1 DID1,   FAIL,   FATS,   FSAVED, HAVDIF, IEND,   INIT,   ROUNDF,
143     2 XCDOBT, PAD
144      SAVE /DINTNC/, /DINTC/
145C
146C     THE VARIABLES HERE DEFINE THE MACHINE ENVIRONMENT.  ALL ARE SET
147C     IN DINTOP.  THE MEANING ATTACHED TO THESE VARIABLES CAN BE
148C     FOUND BY LOOKING AT THE DEFINITIONS IN DINTOP.
149      DOUBLE PRECISION
150     1  EMEPS,  EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF,
151     2  ESMALL, ENZER,  EDELM1, ENINF
152      COMMON /DINTEC/
153     1  EMEPS,  EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF,
154     2  ESMALL, ENZER,  EDELM1, ENINF
155      SAVE /DINTEC/
156C
157C     *****     STATEMENT FUNCTIONS     ********************************
158C
159C SOLVE   PROVIDES THE SOLUTION OF A QUADRATIC EQUATION.
160      SOLVE(SX,SG)=SQRTTB*SX/(SG*SQRTTB+SQRT(ABS(TB)*SG*SG+SX))
161C TDECR   IS USED TO TRANSFORM AN ABSCISSA FROM THE CURRENT COORDINATE
162C         SYSTEM TO ONE IN WHICH NSUB IS DECREMENTED BY A FACTOR OF 2.
163C     TDECR(SX)=TA+(SX-TA)*((SX-TA)/TB)
164      TDECR(SX)=TA*(1.0+TA/TB)+SX*((SX-TA)/TB-TA/TB)
165C
166C     *****     EXECUTABLE STATEMENTS     ******************************
167C
168      SX=SXMIN
169      IF (NSUB .EQ. 0) THEN
170         SG = SX
171      ELSE
172         SG = TDECR(SX)
173         IF (NSUB .NE. 2) SG = TDECR(SG)
174      END IF
175      SMIN=EDELM3*MAX(EDELM1,ABS(SG))
176      IF (NSUB .NE. 0) THEN
177         SQRTTB=SQRT(ABS(TB))
178         SG = ABS((SX-TA)/TB)
179         SMIN=SOLVE(SMIN,SG)
180         IF (NSUB .NE. 2) SMIN=SOLVE(SMIN,SG*SG)
181      END IF
182      DINTSM=SMIN
183      RETURN
184      END
185