**** options ****
jprobl=1 square
2 notch cut out of square
5 squeezed square
6 quarter circle
7 quarter of circular fin waveguide 1,1
11 quarter of circular fin waveguide general
jguess=1 the trivial guess 0. (all cases)
6 an excellent guess for the notch
8 guess based on arg(wt(i)), beta=0
9 assumes x(1),...,x(2*n-4) set in main.f
nrml = 1 0 -> 0, wn -> wtn, side 1 -> circle of side_t 1
3 optimize vertex errors, starting from nrml=1
jeqn = 2 re, im wn = wtn
5 sigma(j) = sigmat(j), 1<=j<=n-1
**** program outline ****
main
sincp
sincp2
scrss
n2f
sfneq
sinteg
sstpar
selim
strans
sray
ode
sfode
smtov
scirc
sle2
snewt2
sinter
scrss
smapc
srinv
srarg
scplot
-----------------------------------------------------------
scirc find radius and orientation from derivatives
scplot plot polygon
scrss search for parameter in normalizing transformation
selim solve for 3 betas
sexprl relative exponential
sfneq define nonlinear system
sfode define ode system
sincp initialize problem and parameters
sincp2 find auxiliary information from initial data
sinteg given parameters, find computed polygon
sinter nonlinear system for intersection of circles
sle2 2 by 2 linear systems
smapc find center of Mobius image of circle
smobs apply implicitly defined Mobius transformation
smtov find vertices from midpoint data
snewt2 newton's method
srarg compute relative angle about center
sray integrate ode to midpoints of sides
srinv zero-tolerant divide
sstpar store new guess at parameters
strans transform to constrained variables
-----------------------------------------------------------
current target
w wt vertices
wc wtc centers
wdir wtdir orientations
side j connects w(j-1) and w(j)
boundary is traced counterclockwise around polygon
wdir = 1 for convex polygons
miscellaneous common variables:
epsfcn estimate of error in computing sfneq
common variables no longer used:
s
punp
ptrace
hybeps
c n number of sides
c wtdir orientations (convex=1, concave=-1)
c w0 image of 0 inside polygon
c wt vertices
c wtm point on side wt(n) -- wt(1)
c wtc centers of sides
c ****** conformal mapping of circular arc polygons ******
c Petter Bjorstad and Eric Grosse 7 Apr 1986
c Siam J Sci Stat Comp Apr 1987?
c
c CAUTION. This is an experimental research code. We certainly
c cannot promise you it will always converge quickly! Please
c let us know what troubles you have; although there are no
c immediate plans to turn it into good mathematical software,
c we would certainly welcome any improvements.
c
c The routine N2F called here is availble by asking netlib:
c send n2f ivset from port
c
real x(16)
c
common /initc/ a, b, c, d, winit1, winit2
complex a, b, c, d, winit1, winit2
common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun
integer maxfun
real radeps,odeeps,hybeps,epsfcn,dmax
c
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd, zk, ci, zero, wc(10), w0
c
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
common /plotcm/ punp, ptrace, mtrace, ntrace, odet
integer punp, ptrace, mtrace, ntrace
complex odet(100,10)
common /method/ nrml, jeqn
integer nrml, jeqn
common /wavegd/ cfws, cfwt
real cfws, cfwt
c
integer m, neqn, nl2iv(2000)
c 82*m
real nl2v(3000)
c 105 + m*(neqn+2*m+17) + 2*neqn
complex z2,wz
external sfode, sfneq
real y(6), r1, r2, relerr, abserr, awork(226)
integer iwork(5),iflag
integer i, j, k, nfun, niter
c
integer mxfcal, mxiter, afctol, rfctol, dltfdj, sctol
integer xctol, delta0, dltfdc, lmax0, covprt, rdreq
parameter ( mxfcal=17, mxiter=18, afctol=31, rfctol=32 )
parameter ( dltfdj=43, sctol=37, xctol=33, delta0=44 )
parameter ( dltfdc=42, lmax0=35, covprt=14, rdreq=57 )
c
cfws=1.1
cfwt=1
c from converged 1,1
x( 1) = -4.412843E-01
x( 2) = -9.410260E-01
x( 3) = 1.542326E+00
x( 4) = -2.231437E+00
x( 5) = 2.489092E-02
x( 6) = -2.438678E-01
call sincp(x,2,1,1,2)
m = 2*n-4
c
if(jeqn.eq.1)then
neqn = 2*n-1
else if(jeqn.eq.2)then
neqn = 2*n-2
else if(jeqn.eq.3)then
neqn = 3*n-2
else if(jeqn.eq.4)then
neqn = n-1
else if(jeqn.eq.5)then
neqn = n-1
else
write(6,1119) jeqn
1119 format(' equation schema ',i5,' is not implemented')
stop
endif
c
it = 0
call ivset(1,nl2iv,2000,3000,nl2v)
nl2iv(mxfcal) = maxfun
nl2iv(mxiter) = 50
c nl2iv(covprt) = 0
c nl2iv(rdreq ) = 0
nl2v(lmax0 ) = dmax
nl2v(afctol) = epsfcn**2
nl2v(rfctol) = epsfcn
nl2v(sctol ) = epsfcn
nl2v(dltfdj) = sqrt(epsfcn)
nl2v(xctol ) = max(0.1*eps**(0.5),10*epsfcn)
nl2v(delta0) = sqrt(epsfcn)
nl2v(dltfdc) = (epsfcn)**(1./3.)
call n2f(neqn,m,x,sfneq,nl2iv,2000,3000,nl2v,nl2iv,nl2v,sfode)
c final function evaluation to be sure common block values are
c up to date
call sfneq(neqn,m,x,i,nl2v,nl2iv,nl2v,sfode)
c write(6,*) 'hessian='
c k=nl2iv(56)-1
c do 569 i = 1, m
c write(6,*) (nl2v(j+k),j=1,i)
c k = k + i
c 569 continue
do 10 i = 1, m
write (6,1003) i, x(i)
1003 format (' x(', i2, ') =', 1pe14.6)
10 continue
abserr = 0
do 15 i = 1, n
abserr = max(abserr,abs(w(i)-wt(i)))
15 continue
write(6,1005) abserr
1005 format('max abs vertex err =',1pe10.2)
do 20 i = 1, n
write (6,1004) i, beta(i)
1004 format ('beta(', i2, ') =', 1pe14.6)
20 continue
call resist(n,x)
c
stop
end
C subroutine scplot(wb, wa, wd)
C integer wd(10)
C complex wb(10), wa(10)
C end
Compliments of netlib Fri Apr 19 13:21:51 EST 1985
c modified to flag finite difference calls to calcr by ifidif=1
SUBROUTINE N2F(N, P, X, CALCR, IV, LIV, LV, V,
1 UIPARM, URPARM, UFPARM)
common /scplt/ ifidif
C
C *** MINIMIZE A NONLINEAR SUM OF SQUARES USING RESIDUAL VALUES ONLY..
C *** THIS AMOUNTS TO N2G WITHOUT THE SUBROUTINE PARAMETER CALCJ.
C
C *** PARAMETERS ***
C
INTEGER N, P, LIV, LV
C/6
INTEGER IV(LIV), UIPARM(1)
REAL X(P), V(LV), URPARM(1)
C/7
C INTEGER IV(LIV), UIPARM(*)
C REAL X(P), V(LV), URPARM(*)
C/
EXTERNAL CALCR, UFPARM
C
C----------------------------- DISCUSSION ----------------------------
C
C THIS AMOUNTS TO SUBROUTINE NL2SNO (REF. 1) MODIFIED TO CALL
C RN2G.
C THE PARAMETERS FOR N2F ARE THE SAME AS THOSE FOR N2G
C (WHICH SEE), EXCEPT THAT CALCJ IS OMITTED. INSTEAD OF CALLING
C CALCJ TO OBTAIN THE JACOBIAN MATRIX OF R AT X, N2F COMPUTES
C AN APPROXIMATION TO IT BY FINITE (FORWARD) DIFFERENCES -- SEE
C V(DLTFDJ) BELOW. N2F USES FUNCTION VALUES ONLY WHEN COMPUT-
C THE COVARIANCE MATRIX (RATHER THAN THE FUNCTIONS AND GRADIENTS
C THAT N2G MAY USE). TO DO SO, N2F SETS IV(COVREQ) TO MINUS
C ITS ABSOLUTE VALUE. THUS V(DELTA0) IS NEVER REFERENCED AND ONLY
C V(DLTFDC) MATTERS -- SEE NL2SOL FOR A DESCRIPTION OF V(DLTFDC).
C THE NUMBER OF EXTRA CALLS ON CALCR USED IN COMPUTING THE JACO-
C BIAN APPROXIMATION ARE NOT INCLUDED IN THE FUNCTION EVALUATION
C COUNT IV(NFCALL), BUT ARE RECORDED IN IV(NGCALL) INSTEAD.
C
C V(DLTFDJ)... V(43) HELPS CHOOSE THE STEP SIZE USED WHEN COMPUTING THE
C FINITE-DIFFERENCE JACOBIAN MATRIX. FOR DIFFERENCES IN-
C VOLVING X(I), THE STEP SIZE FIRST TRIED IS
C V(DLTFDJ) * MAX(ABS(X(I)), 1/D(I)),
C WHERE D IS THE CURRENT SCALE VECTOR (SEE REF. 1). (IF
C THIS STEP IS TOO BIG, I.E., IF CALCR SETS NF TO 0, THEN
C SMALLER STEPS ARE TRIED UNTIL THE STEP SIZE IS SHRUNK BE-
C LOW 1000 * MACHEP, WHERE MACHEP IS THE UNIT ROUNDOFF.
C DEFAULT = MACHEP**0.5.
C
C *** REFERENCE ***
C
C 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE
C NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH.
C SOFTWARE, VOL. 7, NO. 3.
C
C *** GENERAL ***
C
C CODED BY DAVID M. GAY.
C
C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++
C
C *** EXTERNAL SUBROUTINES ***
C
EXTERNAL IVSET, RN2G, N2RDP
C
C IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS.
C RN2G... CARRIES OUT OPTIMIZATION ITERATIONS.
C N2RDP... PRINTS REGRESSION DIAGNOSTICS.
C
C *** LOCAL VARIABLES ***
C
INTEGER D1, DK, DR1, I, IV1, J1K, K, N1, N2, NF, NG, RD1, R1, RN
REAL H, H0, HLIM, NEGPT5, ONE, XK, ZERO
C
C *** IV AND V COMPONENTS ***
C
INTEGER COVREQ, D, DINIT, DLTFDJ, J, MODE, NEXTV, NFCALL, NFGCAL,
1 NGCALL, NGCOV, R, REGD, REGD0, TOOBIG, VNEED
C/6
DATA COVREQ/15/, D/27/, DINIT/38/, DLTFDJ/43/, J/70/, MODE/35/,
1 NEXTV/47/, NFCALL/6/, NFGCAL/7/, NGCALL/30/, NGCOV/53/,
2 R/61/, REGD/67/, REGD0/82/, TOOBIG/2/, VNEED/4/
C/7
C PARAMETER (COVREQ=15, D=27, DINIT=38, DLTFDJ=43, J=70, MODE=35,
C 1 NEXTV=47, NFCALL=6, NFGCAL=7, NGCALL=30, NGCOV=53,
C 2 R=61, REGD=67, REGD0=82, TOOBIG=2, VNEED=4)
C/
DATA HLIM/0.1E+0/, NEGPT5/-0.5E+0/, ONE/1.E+0/, ZERO/0.E+0/
C
C--------------------------------- BODY ------------------------------
C
IF (IV(1) .EQ. 0) CALL IVSET(1, IV, LIV, LV, V)
IV(COVREQ) = -IABS(IV(COVREQ))
IV1 = IV(1)
IF (IV1 .EQ. 14) GO TO 10
IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10
IF (IV1 .EQ. 12) IV(1) = 13
IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+2)
CALL RN2G(X, V, IV, LIV, LV, N, N, N1, N2, P, V, V, V, X)
IF (IV(1) .NE. 14) GO TO 999
C
C *** STORAGE ALLOCATION ***
C
IV(D) = IV(NEXTV)
IV(R) = IV(D) + P
IV(REGD0) = IV(R) + N
IV(J) = IV(REGD0) + N
IV(NEXTV) = IV(J) + N*P
IF (IV1 .EQ. 13) GO TO 999
C
10 D1 = IV(D)
DR1 = IV(J)
R1 = IV(R)
RN = R1 + N - 1
RD1 = IV(REGD0)
C
20 CALL RN2G(V(D1), V(DR1), IV, LIV, LV, N, N, N1, N2, P, V(R1),
1 V(RD1), V, X)
IF (IV(1)-2) 30, 50, 100
C
C *** NEW FUNCTION VALUE (R VALUE) NEEDED ***
C
30 NF = IV(NFCALL)
ifidif = 0
CALL CALCR(N, P, X, NF, V(R1), UIPARM, URPARM, UFPARM)
IF (NF .GT. 0) GO TO 40
IV(TOOBIG) = 1
GO TO 20
40 IF (IV(1) .GT. 0) GO TO 20
C
C *** COMPUTE FINITE-DIFFERENCE APPROXIMATION TO DR = GRAD. OF R ***
C
C *** INITIALIZE D IF NECESSARY ***
C
50 IF (IV(MODE) .LT. 0 .AND. V(DINIT) .EQ. ZERO)
1 CALL V7SCP(P, V(D1), ONE)
C
J1K = DR1
DK = D1
NG = IV(NGCALL) - 1
IF (IV(1) .EQ. (-1)) IV(NGCOV) = IV(NGCOV) - 1
DO 90 K = 1, P
XK = X(K)
H = V(DLTFDJ) * AMAX1( ABS(XK), ONE/V(DK))
H0 = H
DK = DK + 1
60 X(K) = XK + H
NF = IV(NFGCAL)
ifidif = 0
CALL CALCR (N, P, X, NF, V(J1K), UIPARM, URPARM, UFPARM)
NG = NG + 1
IF (NF .GT. 0) GO TO 70
H = NEGPT5 * H
IF ( ABS(H/H0) .GE. HLIM) GO TO 60
IV(TOOBIG) = 1
IV(NGCALL) = NG
GO TO 20
70 X(K) = XK
IV(NGCALL) = NG
DO 80 I = R1, RN
V(J1K) = (V(J1K) - V(I)) / H
J1K = J1K + 1
80 CONTINUE
90 CONTINUE
GO TO 20
C
100 IF (IV(REGD) .GT. 0) IV(REGD) = RD1
CALL N2RDP(IV, LIV, LV, N, V(RD1), V)
C
999 RETURN
C
C *** LAST CARD OF N2F FOLLOWS ***
END
c name changed for use inside sinteg
Compliments of netlib Fri Apr 19 13:21:51 EST 1985
c modified to flag finite difference calls to calcr by ifidif=1
SUBROUTINE N2F2(N, P, X, CALCR, IV, LIV, LV, V,
1 UIPARM, URPARM, UFPARM)
common /scplt/ ifidif
C
C *** MINIMIZE A NONLINEAR SUM OF SQUARES USING RESIDUAL VALUES ONLY..
C *** THIS AMOUNTS TO N2G WITHOUT THE SUBROUTINE PARAMETER CALCJ.
C
C *** PARAMETERS ***
C
INTEGER N, P, LIV, LV
C/6
INTEGER IV(LIV), UIPARM(1)
REAL X(P), V(LV), URPARM(1)
C/7
C INTEGER IV(LIV), UIPARM(*)
C REAL X(P), V(LV), URPARM(*)
C/
EXTERNAL CALCR, UFPARM
C
C----------------------------- DISCUSSION ----------------------------
C
C THIS AMOUNTS TO SUBROUTINE NL2SNO (REF. 1) MODIFIED TO CALL
C RN2G.
C THE PARAMETERS FOR N2F ARE THE SAME AS THOSE FOR N2G
C (WHICH SEE), EXCEPT THAT CALCJ IS OMITTED. INSTEAD OF CALLING
C CALCJ TO OBTAIN THE JACOBIAN MATRIX OF R AT X, N2F COMPUTES
C AN APPROXIMATION TO IT BY FINITE (FORWARD) DIFFERENCES -- SEE
C V(DLTFDJ) BELOW. N2F USES FUNCTION VALUES ONLY WHEN COMPUT-
C THE COVARIANCE MATRIX (RATHER THAN THE FUNCTIONS AND GRADIENTS
C THAT N2G MAY USE). TO DO SO, N2F SETS IV(COVREQ) TO MINUS
C ITS ABSOLUTE VALUE. THUS V(DELTA0) IS NEVER REFERENCED AND ONLY
C V(DLTFDC) MATTERS -- SEE NL2SOL FOR A DESCRIPTION OF V(DLTFDC).
C THE NUMBER OF EXTRA CALLS ON CALCR USED IN COMPUTING THE JACO-
C BIAN APPROXIMATION ARE NOT INCLUDED IN THE FUNCTION EVALUATION
C COUNT IV(NFCALL), BUT ARE RECORDED IN IV(NGCALL) INSTEAD.
C
C V(DLTFDJ)... V(43) HELPS CHOOSE THE STEP SIZE USED WHEN COMPUTING THE
C FINITE-DIFFERENCE JACOBIAN MATRIX. FOR DIFFERENCES IN-
C VOLVING X(I), THE STEP SIZE FIRST TRIED IS
C V(DLTFDJ) * MAX(ABS(X(I)), 1/D(I)),
C WHERE D IS THE CURRENT SCALE VECTOR (SEE REF. 1). (IF
C THIS STEP IS TOO BIG, I.E., IF CALCR SETS NF TO 0, THEN
C SMALLER STEPS ARE TRIED UNTIL THE STEP SIZE IS SHRUNK BE-
C LOW 1000 * MACHEP, WHERE MACHEP IS THE UNIT ROUNDOFF.
C DEFAULT = MACHEP**0.5.
C
C *** REFERENCE ***
C
C 1. DENNIS, J.E., GAY, D.M., AND WELSCH, R.E. (1981), AN ADAPTIVE
C NONLINEAR LEAST-SQUARES ALGORITHM, ACM TRANS. MATH.
C SOFTWARE, VOL. 7, NO. 3.
C
C *** GENERAL ***
C
C CODED BY DAVID M. GAY.
C
C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++
C
C *** EXTERNAL SUBROUTINES ***
C
EXTERNAL IVSET, RN2G, N2RDP
C
C IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS.
C RN2G... CARRIES OUT OPTIMIZATION ITERATIONS.
C N2RDP... PRINTS REGRESSION DIAGNOSTICS.
C
C *** LOCAL VARIABLES ***
C
INTEGER D1, DK, DR1, I, IV1, J1K, K, N1, N2, NF, NG, RD1, R1, RN
REAL H, H0, HLIM, NEGPT5, ONE, XK, ZERO
C
C *** IV AND V COMPONENTS ***
C
INTEGER COVREQ, D, DINIT, DLTFDJ, J, MODE, NEXTV, NFCALL, NFGCAL,
1 NGCALL, NGCOV, R, REGD, REGD0, TOOBIG, VNEED
C/6
DATA COVREQ/15/, D/27/, DINIT/38/, DLTFDJ/43/, J/70/, MODE/35/,
1 NEXTV/47/, NFCALL/6/, NFGCAL/7/, NGCALL/30/, NGCOV/53/,
2 R/61/, REGD/67/, REGD0/82/, TOOBIG/2/, VNEED/4/
C/7
C PARAMETER (COVREQ=15, D=27, DINIT=38, DLTFDJ=43, J=70, MODE=35,
C 1 NEXTV=47, NFCALL=6, NFGCAL=7, NGCALL=30, NGCOV=53,
C 2 R=61, REGD=67, REGD0=82, TOOBIG=2, VNEED=4)
C/
DATA HLIM/0.1E+0/, NEGPT5/-0.5E+0/, ONE/1.E+0/, ZERO/0.E+0/
C
C--------------------------------- BODY ------------------------------
C
write(6,*) 'n2f2{'
IF (IV(1) .EQ. 0) CALL IVSET(1, IV, LIV, LV, V)
IV(COVREQ) = -IABS(IV(COVREQ))
IV1 = IV(1)
IF (IV1 .EQ. 14) GO TO 10
IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10
IF (IV1 .EQ. 12) IV(1) = 13
IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+2)
CALL RN2G(X, V, IV, LIV, LV, N, N, N1, N2, P, V, V, V, X)
IF (IV(1) .NE. 14) GO TO 999
C
C *** STORAGE ALLOCATION ***
C
IV(D) = IV(NEXTV)
IV(R) = IV(D) + P
IV(REGD0) = IV(R) + N
IV(J) = IV(REGD0) + N
IV(NEXTV) = IV(J) + N*P
IF (IV1 .EQ. 13) GO TO 999
C
10 D1 = IV(D)
DR1 = IV(J)
R1 = IV(R)
RN = R1 + N - 1
RD1 = IV(REGD0)
C
20 CALL RN2G(V(D1), V(DR1), IV, LIV, LV, N, N, N1, N2, P, V(R1),
1 V(RD1), V, X)
IF (IV(1)-2) 30, 50, 100
C
C *** NEW FUNCTION VALUE (R VALUE) NEEDED ***
C
30 NF = IV(NFCALL)
ifidif = 0
CALL CALCR(N, P, X, NF, V(R1), UIPARM, URPARM, UFPARM)
IF (NF .GT. 0) GO TO 40
IV(TOOBIG) = 1
GO TO 20
40 IF (IV(1) .GT. 0) GO TO 20
C
C *** COMPUTE FINITE-DIFFERENCE APPROXIMATION TO DR = GRAD. OF R ***
C
C *** INITIALIZE D IF NECESSARY ***
C
50 IF (IV(MODE) .LT. 0 .AND. V(DINIT) .EQ. ZERO)
1 CALL V7SCP(P, V(D1), ONE)
C
J1K = DR1
DK = D1
NG = IV(NGCALL) - 1
IF (IV(1) .EQ. (-1)) IV(NGCOV) = IV(NGCOV) - 1
DO 90 K = 1, P
XK = X(K)
H = V(DLTFDJ) * AMAX1( ABS(XK), ONE/V(DK))
H0 = H
DK = DK + 1
60 X(K) = XK + H
NF = IV(NFGCAL)
ifidif = 0
CALL CALCR (N, P, X, NF, V(J1K), UIPARM, URPARM, UFPARM)
NG = NG + 1
IF (NF .GT. 0) GO TO 70
H = NEGPT5 * H
IF ( ABS(H/H0) .GE. HLIM) GO TO 60
IV(TOOBIG) = 1
IV(NGCALL) = NG
GO TO 20
70 X(K) = XK
IV(NGCALL) = NG
DO 80 I = R1, RN
V(J1K) = (V(J1K) - V(I)) / H
J1K = J1K + 1
80 CONTINUE
90 CONTINUE
GO TO 20
C
100 IF (IV(REGD) .GT. 0) IV(REGD) = RD1
CALL N2RDP(IV, LIV, LV, N, V(RD1), V)
C
999 continue
write(6,*) 'n2f2}'
return
END
subroutine resist ( n, x )
integer n
real x(25), ang(25), ang2(25)
complex z(25)
integer i,nm
real k, r
complex a,b,c,d,chi
c
call stransa(x,z,ang,n)
write(6,1001)(i,ang(i),i=1,n)
1001 format('tau(',i2,')=',f10.6)
write(6,1002)(i,z(i),i=1,n)
1002 format('z(',i2,')=',1p2e16.8)
c
do 190 i=1,n-3
do 180 j=i+2,n-1
a = z(i)
b = z(i+1)
c = z(j)
d = z(j+1)
chi = (b-a)*(d-c)/((c-b)*(a-d))
k = real(chi)
k = 1 - 2*(k+sqrt(k*k-k))
r = 2*ellip(k)/ellip(sqrt(1-k*k))
write(6,1004) i, j, r
180 continue
190 continue
1004 format('resistance from side',i2,' to side ',i2,' = ',1pe16.8)
return
end
c-------------------------------------------------
real function ellip ( k )
c complete elliptic integral
c from Hart, Computer Approximations, 154
real eps, k, a, b, t, pi
pi = 4*atan(1.)
eps = (r1mach(4))
a = 1
b = sqrt(1-k*k)
10 continue
t = (a+b)/2
b = sqrt(a*b)
a = t
if(abs(a-b).gt.eps*(a+b)) goto 10
ellip = pi/(2*a)
return
end
function sarc ( z1, z2, c, wdir )
c
c arc length between points z1 and z2 on circle centered at c
c wdir = 1 counterclockwise from z1 to z2
c -1 clockwise
c
real sarc
complex z1, z2, c, c1, t
integer wdir
real r, s, pi
pi = 4.*atan(1.)
r = .5*(abs(z1-c)+abs(z2-c))
s = abs(z1-z2)
c phi = angle counterclockwise from z1-c to z2-c
t = z1 - c
t = (z2-c)/t
phi = atan2(aimag(t),real(t))
if(phi.lt.0) phi = phi + 2*pi
if(phi.eq.0.0)then
t = c-z1
c1 = .5*(z1+z2) + 30*max(abs(z1),abs(z2))*t/abs(t)
t = z1 - c1
t = (z2-c1)/t
phi = atan2(aimag(t),real(t))
if(phi.lt.0) phi = phi + 2*pi
endif
s = min(1.,max(-1.,.5*s/r))
if( (wdir.eq. 1 .and. phi.le.pi) .or.
+ (wdir.eq.-1 .and. phi.gt.pi) )then
sarc = 2 * r * asin(s)
else
sarc = 2 * r * (pi-asin(s))
end if
return
end
function oarc ( z1, z2, c)
c
c arc length between points z1 and z2 on circle centered at c
c defined in some magic way ?!?!
c
real oarc
complex z1, z2, c
real pi, arg1, arg2, t, t0
pi=4*atan(1.e0)
arg1= atan2(aimag(z1-c), real(z1-c))
arg2= atan2(aimag(z2-c), real(z2-c))
t=arg1-arg2
if(abs(z1-z2).le.0.1e0*abs(z2-c)) t=asin(aimag((z1-z2)/(z2-c)))
t=mod(t+3*pi,2*pi)-pi
if(abs(t).lt.0.01)goto 100
oarc=abs(z2-c)*t
return
100 t0=t
t=t*t
oarc=abs(z1-z2)/sqrt(1-(t/(3*4))*(1-(t/(5*6))*(1-(t/(7*8))*(1-
$ (t/(9*10))*(1-(t/(11*12))*(1-(t/(13*14))))))))
if(t0.lt.0) oarc = - oarc
return
end
subroutine scirc(z, w, w1, w2, c, r, orient)
c
integer orient
real r
complex z, w, w1, w2, c
c
c this routine computes the radius, and the
c orientation of the arc from derivative information.
c
real phidot, w11r, w11i, w11s, w22r, w22i
complex z1, w11, w22
c
z1 = z* cmplx(0.e0, 1.e0)
w11 = w1*z1
w11r = real(w11)
w11i = aimag(w11)
w11s = w11r**2+w11i**2
w22 = -z*(z*w2+w1)
w22r = real(w22)
w22i = aimag(w22)
phidot = w11r*w22i-w11i*w22r
orient = sign(1.e0,phidot)
r = abs(phidot)/w11s
r = 1.e0/r
c = w + cmplx(0.e0,1.e0)*orient*r*w11
r = r * sqrt(w11s)
return
end
subroutine scplot(wb, wa, wd)
c
c trace boundary of circular arc polygon
c
integer wd(10)
complex wb(10), wa(10)
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd, zk, ci, zero, wc(10), w0
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
common /scplt/ ifidif
integer i
c
if(ifidif.eq.1) return
do 10 i = 1, n
write (6,1002) wb(i), wa(i), wd(i)
1002 format ('* ',1p4e15.7,i3)
10 continue
return
end
subroutine scrss(wm,w1,wn,aa)
c
complex wm,w1,wn,aa
c
c this routine finds aa , the parameter in an intermediate
c mobius transformation that enforces the circle to
c circle condition.
c
common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun
integer maxfun
real radeps,odeeps,hybeps,epsfcn,dmax
c
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
c
real a,b,c,d,dmb,cma,cl,dnom
real x,y
c
a = real(wm/wn)
b = aimag(wm/wn)
c = real(w1/wn)
d = aimag(w1/wn)
c
dnom = (a-1.e0)*(a+1.e0) + b*b
if( abs(dnom).lt.radeps) go to 10
a = (a-1.e0)/dnom
b = b/dnom
dnom = (c-1.e0)*(c+1.e0) + d*d
if( abs(dnom).lt.radeps) go to 20
c = (c-1.e0)/dnom
d = d/dnom
c
dmb = d-b
cma = c-a
dnom = dmb*dmb+cma*cma
x = (cma*(cma-2.e0*b*dmb)+dmb*dmb*(2.e0*a-1.e0))/dnom
y = 2.e0*cma*(b*cma+(1.e0-a)*dmb)/dnom
aa= cmplx(x,y)
return
10 aa= cmplx(1.e0,0.e0)
cl = (a-1.e0)/b
dnom = (c-1.e0)*(c+1.e0) + d*d
if( abs(dnom).lt.radeps) return
c = (c-1.e0)/dnom
d = d/dnom
x = (cl*(cl-2.e0*d)+2.e0*c-1.e0)/(cl*cl+1.e0)
y = 2.e0*cl*(d*cl-c+1.e0)/(cl*cl+1.e0)
aa= cmplx(x,y)
return
20 cl = (c-1.e0)/d
x = (cl*(cl-2.e0*b)+2.e0*a-1.e0)/(cl*cl+1.e0)
y = 2.e0*cl*(b*cl-a+1.e0)/(cl*cl+1.e0)
aa= cmplx(x,y)
return
end
subroutine selim(n,a,bet,gam)
c
integer n,i
real a(1),bet(1),gam(1)
c
c eliminates algebraic sideconditions.
c
c input.
c a(i) i=1,n -vertex coordinates in halfplane.
c (on real axis.)
c bet(i) i=1,n-3 -first n-3 beta-values.
c gam(i) i=1,n -given constants.
c output.
c bet(n-2),bet(n-1) and bet(n).
c
real a1,a2,a3,b1,b2,b3
c
b1=0.e0
b2=0.e0
b3=0.e0
bet(n-2)=0.e0
bet(n-1)=0.e0
bet(n) =0.e0
c
c compute right hand side.
c
do 10 i=1,n
b1=b1+bet(i)
b2=b2+2.e0*bet(i)*a(i)+gam(i)
b3=b3+a(i)*(bet(i)*a(i)+gam(i))
10 continue
a1=a(n-2)
a2=a(n-1)
a3=a(n)
b1=b1-(bet(n)+bet(n-1)+bet(n-2))
b2=b2-2.e0*(bet(n)*a3+bet(n-1)*a2+bet(n-2)*a1)
b3=b3-(bet(n)*a3*a3+bet(n-1)*a2*a2+bet(n-2)*a1*a1)
c
c start elimination.
c
b1=-b1
b2=-(b2+2.e0*b1*a1)
b3=-(b3+.5e0*b2*(a1+a2)+b1*a1*a1)
b3=b3/((a3-a1)*(a3-a2))
b2=.5e0*(b2-2.e0*b3*(a3-a1))/(a2-a1)
b1=b1-b2-b3
bet(n-2)=b1
bet(n-1)=b2
bet(n) =b3
return
end
function sexprl (z)
c complex version of Wayne's Aug 80 fnlib routine eric 3 Apr 81
c simplified error handling
c undeclare abs, cexp -> exp 18 apr 85
c
c evaluate (cexp(z)-1)/z . for small abs(z), we use the taylor
c series. we could instead use the expression
c sexprl(z) = (exp(x)*exp(i*y)-1)/z
c = (x*exprel(x) * (1 - 2* sin(y/2)**2) - 2* sin(y/2)**2
c + i* sin(y)*(1+x*exprel(x))) / z
c
complex z, sexprl
real rbnd, sqeps, drlz, r1mach
real r, xn, xln, alneps
data nterms, rbnd, sqeps / 0, 2*0.0e0 /
c
if (nterms.ne.0) go to 10
alneps = log(r1mach(3))
xn = 3.72 - 0.3*alneps
xln = log((xn+1.0)/1.36)
nterms = xn - (xn*xln+alneps)/(xln+1.36) + 1.5
rbnd = r1mach(3)
sqeps = 0.0
c
10 r = abs(z)
if (r.le.0.5) go to 20
c
sexprl = (exp(z) - 1.0) / z
drlz= real(z)
if( abs(drlz).lt.sqeps.and.abs(sexprl).lt.sqeps) goto 101
return
c
20 sexprl = (1.0, 0.0)
if (r.lt.rbnd) return
c
sexprl = (0.0, 0.0)
do 30 i=1,nterms
sexprl = 1.0 + sexprl*z/ float(nterms+2-i)
30 continue
return
c
101 write(6,1001)
1001 format('zexprl answer lt half prec, z near ni2pi')
stop
end
subroutine sfneq(neqn, m, x, nf, f, nl2iv, nl2v, foo)
c
integer m, neqn, nl2iv(1)
real x(m), f(neqn), nl2v(1)
external foo
c
c this routine defines the system of nonlinear equations.
c
common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun
integer maxfun
real radeps,odeeps,hybeps,epsfcn,dmax
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd, zk, ci, zero, wc(10), w0
c
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
common /method/ nrml, jeqn
integer nrml, jeqn
c
common /sfail/ ifail
integer ifail
c
integer i, k, is
it=it+1
ifail=0
call sinteg(x)
if(ifail.eq.0)goto 100
nf=0
return
100 continue
if(jeqn.eq.1)then
do 111 i = 1, n-1
k = 2*i
f(k-1) = real(wt(i)-w(i))
f(k) = aimag(wt(i)-w(i))
111 continue
do 112 i=1,1
f(2*n-2+i) = arct(i) - arcc(i)
112 continue
else if(jeqn.eq.2)then
do 121 i = 1, n-1
k = 2*i
f(k-1) = real(wt(i)-w(i))
f(k) = aimag(wt(i)-w(i))
121 continue
else if(jeqn.eq.3)then
do 131 i = 1, n-1
k = 2*i
f(k-1) = real(wt(i)-w(i))
f(k) = aimag(wt(i)-w(i))
131 continue
do 132 i=1,n
f(2*n-2+i) = arct(i) - arcc(i)
132 continue
else if(jeqn.eq.4)then
do 141 i = 1, n-1
f(i) = abs(wt(i)-w(i))
141 continue
else if(jeqn.eq.5)then
do 151 i = 1, n-1
f(i) = arct(i) - arcc(i)
151 continue
end if
call scplot(w,wc,wdir)
return
end
subroutine sfode(r,y,yp)
c
real r,y(6),yp(6)
c
c routine called by ode. direct calculation.
c for the schwartzian , third order equation.
c
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd,zk,ci,zero, wc(10), w0
c
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
c
common/node/ nfun
integer nfun
c
integer i
real ar, ai, br, bi
real denom
real z2r, z2i, csumr, csumi, c2r, c2i, zkr, zki
real cr, cii, zdr, zdi, c1r, c1i
c
c zk is cut in circle-halfplane map.
c zd is exp(%i*teta), teta direction of integration.
c z(i) , i=1,2, n. contains the vertices in the
c unit disk. (singular points.)
c
c--------------
c for debugging
if(r.lt.0)then
write(6,*)" zk=cmplx",zk
write(6,*)" zd=cmplx",zd
write(6,*)" z2=",z2
write(6,*)" ci=cmplx",ci
write(6,*)" c2=",c2
do 8881 i=1,n
write(6,*)" cc(",i,")=cmplx",cc(i)
write(6,*)" gamma(",i,")=",gamma(i)
write(6,*)" beta(",i,")=",beta(i)
8881 continue
return
end if
c--------------
nfun = nfun + 1
zkr = real(zk)
zki = aimag(zk)
zdr = real(zd)
zdi = aimag(zd)
c cy2= cmplx(y(3),y(4))
c cy3= cmplx(y(5),y(6))
z2r=r*zdr
z2i=r*zdi
c csum = zero
csumr = 0.e0
csumi = 0.e0
c c2 = 1.e0/(zk-z2)
denom = (zkr-z2r)**2 + (zki-z2i)**2
c2r = (zkr-z2r) / denom
c2i = - (zki-z2i) / denom
c c = ci*c2*(zk+z2)
ar = zkr+z2r
ai = zki+z2i
cii = c2r*ar-c2i*ai
cr = -(c2i*ar+c2r*ai)
c c2 = 8.e0*(c2*c2*zd*zk)**2
ar = c2r*c2r-c2i*c2i
ai = 2*c2r*c2i
br = ar*zdr-ai*zdi
bi = ai*zdr+ar*zdi
ar = br*zkr-bi*zki
ai = bi*zkr+br*zki
c2r = 8*(ar*ar-ai*ai)
c2i = 8*(2*ar*ai)
do 10 i = 1,n
c c1 = .5e0/(c-cc(i))
ar = cr- real(cc(i))
denom = ar**2+cii**2
c1r = .5e0*ar/denom
c1i = -.5e0*cii/denom
c csum = csum + c1*(c1*gamma(i)+beta(i))
ar = c1r*gamma(i)+beta(i)
ai = c1i*gamma(i)
csumr = csumr + c1r*ar-c1i*ai
csumi = csumi + c1i*ar+c1r*ai
10 continue
c c1 = 1.5e0*cy3**2/cy2-c2*cy2*csum
ar = 1.5e0*(y(5)*y(5)-y(6)*y(6))
ai = 1.5e0*2*y(5)*y(6)
denom = y(3)**2+y(4)**2
c1r = (ar*y(3)+ai*y(4))/denom
c1i = (ai*y(3)-ar*y(4))/denom
ar = c2r*y(3)-c2i*y(4)
ai = c2i*y(3)+c2r*y(4)
c1r = c1r - (ar*csumr-ai*csumi)
c1i = c1i - (ai*csumr+ar*csumi)
yp(1)=y(3)
yp(2)=y(4)
yp(3)=y(5)
yp(4)=y(6)
yp(5)=c1r
yp(6)=c1i
c write(6,*)'r,yp6=',r,yp(6)
return
end
subroutine sincp ( x, jprobl, jguess, inrml, ieqn )
c
c initialize parameters describing target circular arc polygon
c
real x(1)
integer jprobl, jguess
common /initc/ a, b, c, d, winit1, winit2
complex a, b, c, d, winit1, winit2
common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun
integer maxfun
real radeps,odeeps,hybeps,epsfcn,dmax
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd, zk, ci, zero, wc(10), w0
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
common /plotcm/ punp, ptrace, mtrace, ntrace, odet
integer punp, ptrace, mtrace, ntrace
complex odet(100,10)
common /method/ nrml, jeqn
integer nrml, jeqn
integer i, j, k, m
real cfwa, cfws, cfwt,t,t1, h
real th(11)
complex wtm, u
common /wavegd/ cfws, cfwt
c
pi = 4.e0*atan(1.e0)
eps = r1mach(4)
iv(1) = 4
ntrace = 2
nrml = inrml
jeqn = ieqn
c ---- tolerances ----
radeps = 100*eps
odeeps = 2000*eps
epsfcn = 10*odeeps
maxfun = 800
dmax = 1.e-2
hybeps = 10*odeeps
c
ci = cmplx(0.e0, 1.e0)
zero = cmplx(0.e0, 0.e0)
zk = cmplx( cos(0.33333e0), sin(0.33333e0))
c
c ----------- particular regions ----------------
c
c n number of sides
c wtdir orientations (convex=1, concave=-1)
c w0 image of 0 inside polygon
c wt vertices
c wtm point on side wt(n) -- wt(1)
c wtc centers of sides
c
go to (10,20,30,40,50,60,70,30,30,30,710), jprobl
c
10 n = 4
write(6,2010)
2010 format(' square example')
do 15 i = 1, n
wtdir(i) = 1
15 continue
w0 = zero
wt(1) = cmplx( .5e0, .5e0)
wt(2) = cmplx(-.5e0, .5e0)
wt(3) = cmplx(-.5e0, -.5e0)
wt(4) = cmplx( .5e0, -.5e0)
wtc(1) = cmplx(-1.e30, 0.)
wtc(2) = cmplx( 0., -1.e30)
wtc(3) = cmplx( 1.e30, 0.)
wtc(4) = cmplx( 0., 1.e30)
wtm = cmplx( .5e0, .0e0)
c the following is needed to avoid "stiffness" in example 1,1,1,5
zk = cmplx( cos(0.11111e0), sin(0.11111e0))
go to 99
c
20 n = 5
write(6,2020)
2020 format(' notch example')
do 25 i = 1, n
wtdir(i) = 1
25 continue
wtdir(1)=-1
w0 = zero
wt(1) = cmplx( .0e0, .5e0)
wt(2) = cmplx(-.5e0, .5e0)
wt(3) = cmplx(-.5e0, -.5e0)
wt(4) = cmplx( .5e0, -.5e0)
wt(5) = cmplx( .5e0, .0e0)
wtc(1) = cmplx( .5e0, .5e0)
wtc(2) = cmplx( 0., -1.e30)
wtc(3) = cmplx( 1.e30, 0.)
wtc(4) = cmplx( 0., 1.e30)
wtc(5) = cmplx(-1.e30, 0.)
t = .25e0 * (2.e0-sqrt(2.e0))
wtm = cmplx(t,t)
go to 99
c
c
30 stop
40 stop
50 n = 4
write(6,2060)
2060 format(' sqeezed square example')
wtdir(1) = -1
wtdir(2) = 1
wtdir(3) = -1
wtdir(4) = 1
w0 = zero
wt(1) = cmplx( .5e0, .5e0)
wt(2) = cmplx(-.5e0, .5e0)
wt(3) = cmplx(-.5e0, -.5e0)
wt(4) = cmplx( .5e0, -.5e0)
wtc(1) = cmplx(4., 0.)
wtc(2) = cmplx(0., -4.)
wtc(3) = cmplx(-4., 0.)
wtc(4) = cmplx(0., 4.)
wtm = cmplx( real(wtc(1))-abs(wt(1)-wtc(1)), 0.e0)
go to 99
60 n = 3
write(6,2050)
2050 format(' quarter circle example')
wtdir(1) = 1
wtdir(2) = 1
wtdir(3) = 1
w0 = cmplx(.4e0,.4e0)
wt(1) = cmplx( 0.e0, 1.e0) - w0
wt(2) = cmplx( 0.e0, 0.e0) - w0
wt(3) = cmplx( 1.e0, 0.e0) - w0
wtc(1) = cmplx(0., 0.) - w0
wtc(2) = cmplx(1.e30, 0.) - w0
wtc(3) = cmplx(0., 1.e30) - w0
wtm = cmplx( 1/sqrt(2.e0), 1/sqrt(2.e0) ) - w0
w0 = 0
go to 99
c
70 n = 5
cfws = 1
cfwt = 1
c jump here from 710
71 continue
do 72 i = 1, n
wtdir(i) = 1
72 continue
write(6,2070) cfws, cfwt
2070 format(' quarter circular fin waveguide example'
+ /' 1/s, 1/t = ',2f15.2)
cfwa = 1.e0
cfws = 1.e0/cfws
cfwt = 1.e0/cfwt
w0 = cmplx( cfwt/200+cfwa/2, cfws/4 )
wt(1)= cmplx( cfwt/2, cfwa ) - w0
wt(2)= cmplx( cfwt/2, cfws/2 ) - w0
wt(3)= cmplx( 0.e0, cfws/2 ) - w0
wt(4)= 0 - w0
wt(5)= cmplx( cfwt/2+cfwa, 0. ) - w0
wtc(1)=cmplx(cfwt/2,0.) - w0
wtm=cmplx( cfwt/2, 0.e0 ) + cfwa*cmplx(1/sqrt(2.),1/sqrt(2.)) - w0
w0 = 0
do 75 i = 1, 4
u = wt(i+1) - wt(i)
u = u / abs(u)
u = wt(i) - u*(real(u)*real(wt(i))+aimag(u)*aimag(wt(i)))
wtc(i+1) = -1.e10 * u
75 continue
go to 99
710 n = 5
go to 71
c 1 2 3 4 5 6 7 8 9
99 go to (100,110,120,130,140,150,160,170,200),jguess
100 continue
m = 2*n-4
do 105 i=1,m
x(i) = 0.e0
105 continue
go to 200
110 stop
x(1) = 2.3
x(2) = -1.3
x(3) = 0.
x(4) = 1.3
go to 200
120 stop
130 stop
140 stop
150 continue
x( 1) = 2.343616E+00
x( 2) = -1.331171E+00
x( 3) = -6.384445E-04
x( 4) = 1.332062E+00
x( 5) = 1.343304E+00
x( 6) = -6.599826E-01
go to 200
160 stop
170 continue
m = 2*n-4
do 171 i=1,m
x(i) = 0.e0
171 continue
m = n-1
th(1) = 0
th(n+1) = 2*pi
do 173 i=1,m
u = (wt(i)-w0)/(wt(n)-w0)
th(i+1) = atan2(aimag(u),real(u))
if(th(i+1).lt.0) th(i+1) = th(i+1) + 2*pi
173 continue
do 175 i=1,m
x(i) = log( (th(i+1)-th(i))/(th(i+2)-th(i+1)) )
175 continue
goto 200
c
200 call sincp2(x,wtm)
write(6,1002) n
1002 format('* ',i3)
call scplot(wt,wtc,wtdir)
do 210 i = 1, n
write(6,*) 'wt,wtc=', wt(i), wtc(i)
210 continue
write(6,1001) radeps,odeeps
1001 format(
x /' radeps = ',1pe11.3
x /' odeeps = ',1pe11.3
x )
return
end
subroutine sincp2 ( x, wtm )
c
real x(1)
complex wtm
common /initc/ a, b, c, d, winit1, winit2
complex a, b, c, d, winit1, winit2
common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun
integer maxfun
real radeps,odeeps,hybeps,epsfcn,dmax
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd, zk, ci, zero, wc(10), w0
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
common /method/ nrml, jeqn
integer nrml, jeqn
integer i
real sarc, phij, phii
complex bear
c
a = cmplx(1.e0,0.e0)
b = zero
c = zero
d = a
winit1 = a
winit2 = 2.e0*(1.e0/zk-ci)
c
c compute normalizing constant for true polygon
c (this is the last time wtm is used)
if(nrml.eq.1.or.nrml.eq.3)then
call scrss(wtm,wt(1),wt(n),cnr)
else if(nrml.eq.2)then
else
stop
end if
do 350 i = 1, n
j = mod(i,n)+1
bear = wtc(i)-wt(i)
bear = bear / abs(bear)
phii = atan2(aimag(bear),real(bear))
if(wtdir(i).ne.1) phii = phii + pi
bear = wtc(j)-wt(i)
bear = bear / abs(bear)
phij = atan2(aimag(bear),real(bear))
if(wtdir(j).ne.1) phij = phij + pi
alpha(i) = mod(5*pi+phii-phij,2*pi)
gamma(i) = 1.e0-(alpha(i)/pi)**2
350 continue
do 360 i = 1, n
j = mod(i-2+n,n)+1
radt(i) = abs(wt(i)-wtc(i))
arct(i) = sarc(wt(j),wt(i),wtc(i),wtdir(i))
360 continue
return
end
subroutine sinteg(x)
c
real x(1)
c
c this routine is the main driver defining the nonlinear
c system of equations.
c
common /initc/ a, b, c, d, winit1, winit2
complex a, b, c, d, winit1, winit2
c
common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun
integer maxfun
real radeps,odeeps,hybeps,epsfcn,dmax
c
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd, zk, ci, zero, wc(10), w0
c
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
c
common /capcm3/ wm
complex wm(10)
common/node/ nfun
integer nfun
common /plotcm/ punp, ptrace, mtrace, ntrace, odet
integer punp, ptrace, mtrace, ntrace
complex odet(100,10)
common /method/ nrml, jeqn
integer nrml, jeqn
c
common /sfail/ ifail
integer ifail
c
c
integer i, is, j
real marg(10), zarg(10)
real sarc
complex mc, wd(2, 10)
complex aa, mapwn
c
c ...for nrml=3 optimization of normalization by n2f2...
integer nl2iv(300)
real nl2v(300)
external werr
integer rdreq, lmax0
complex copt, dopt
parameter ( rdreq=57, lmax0=35 )
c
c
c write(6,1010) it
1010 format(/' -------------------------- ',i5)
nfun = 0
call sstpar(x)
do 10 i = 1, n
zarg(i)=atan2(aimag(z(i)), real(z(i)))
if(zarg(i).lt.0) zarg(i) = zarg(i) + 2*pi
10 continue
marg(1) = zarg(1)/2.e0
is = n-1
do 20 i = 2, is
marg(i) = zarg(i-1)+(zarg(i)-zarg(i-1))/2.e0
20 continue
marg(n) = zarg(n-1)+(2.e0*pi-zarg(n-1))/2.e0
winit2 = (a/d**2)*(d*winit2-2.e0*c*winit1**2)
winit1 = (a/d)*winit1
c write(6,1019) winit1, winit2
1019 format(' (sinteg) winit1 = ',2e15.7
+ /' winit2 = ',2e15.7 )
mtrace = n
do 30 i = 1, n
call sray(marg(i),wm(i),wd(1,i),wd(2,i),
+ ntrace,odet(1,i))
if(ifail.ne.0) return
30 continue
c write(6,2) nfun
2 format(1x,'ode used ',i7,2x,'function evaluations')
call smtov(wm, wd)
if(ifail.ne.0) return
c
c **** normalization ****
c....
if(nrml.eq.1.or.nrml.eq.3)then
c a0 = 1./w(n)
c a1 = aa
c a2 = cnr
c a3 = 1./wt(n)
call scrss(wm(1),w(1),w(n),aa)
if( (abs(aa-1).le.radeps) .and. (abs(cnr-1).le.radeps) ) then
a = 1./w(n)
b = zero
c = zero
d = 1./wt(n)
else if( (abs(aa-1).le.radeps) )then
a = (1.-aa)/w(n)
b = zero
c = 1./(w(n)*wt(n))
d = aa/wt(n)
else if( (abs(cnr-1).le.radeps) )then
a = aa/w(n)
b = zero
c = 1./(w(n)*wt(n))
d = (aa-1.)/wt(n)
else
a = aa*(1.-cnr)/w(n)
b = zero
c = (aa-cnr)/(w(n)*wt(n))
d = (cnr/wt(n))*(1.-aa)
end if
c....
else if(nrml.eq.2)then
a = 1.
b = zero
aa = wt(1)*wt(n)*(w(n)-w(1))
c = ( wt(1)*w(n) - w(1)*wt(n) ) / aa
d = w(1)*w(n)*(wt(n)-wt(1)) / aa
c....
else
write(6,1038) nrml
1038 format(' normalization ',i5,' is not implemented')
stop
end if
c....
if(nrml.eq.3)then
copt = 0
call ivset(1,nl2iv,300,300,nl2v)
nl2iv(rdreq) = 0
nl2v(lmax0) = .2
call n2f2(2*(n-1),2,copt,werr,
$ nl2iv,300,300,nl2v,nl2iv,nl2v,werr)
dopt = 1 - copt*w(n)
c = copt*a+dopt*c
d = dopt*d
end if
c....
mapwn = a*w(n)/(c*w(n)+d)
c write(6,1040) wm(1),w(1),w(n),aa,cnr,mapwn
1040 format(' wm =',2e13.4
+ /' w1 =',2e13.4
+ /' wn =',2e13.4
+ /' aa =',2e13.4
+ /' cnr =',2e13.4
+ /' mapwn =',2e13.4)
c write(6,1039) a, c, d
1039 format(' a =',2e13.4
+ /' c =',2e13.4
+ /' d =',2e13.4)
c
c apply the mobius transformation
c
c write(6,1041) (w(i),wc(i),i=1,n)
c do 44 i = 1, mtrace
c do 43 j = 1, ntrace
c odet(j,i) = a*odet(j,i)/(c*odet(j,i)+d)
c 43 continue
c 44 continue
do 45 i = 1, n
w(i) = a*w(i)/(c*w(i)+d)
mc = a*wc(i)/(c*wc(i)+d)
call smapc(wc(i), abs(wm(i)-wc(i)), a, b, c, d, wc(i))
wm(i) = a*wm(i)/(c*wm(i)+d)
if (abs(wm(i)-wc(i)) .lt. abs(mc-wc(i))) wdir(i) = -wdir(i)
45 continue
c write(6,1041) (w(i),wc(i),i=1,n)
1041 format(' w, wc =',10(/1x,4e13.4))
c
c restore initialization values to naive guess,
c for debugging
c
c a = cmplx(1.e0,0.e0)
c b = zero
c c = zero
c d = a
c winit1 = a
c winit2 = 2.e0*(1.e0/zk-ci)
c
c compute radius of final circular arcs.
c
do 50 i = 1, n
radc(i) = abs(wc(i)-w(i))
50 continue
c
c compute arc-length of final circular arcs.
c
do 60 i = 1, n
is = mod(i-2+n,n)+1
arcc(i) = sarc(w(is),w(i),wc(i),wdir(i))
60 continue
return
end
c***************************************************************
subroutine werr ( nn, mm, copt, nf, f, nl2iv, nl2v, foo )
common /initc/ a, b, c, d, winit1, winit2
complex a, b, c, d, winit1, winit2
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd, zk, ci, zero, wc(10), w0
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
c
complex copt, dopt, ws, lc, ld
integer i, k, nn, mm, nf, nl2iv(1)
real f(nn), nl2v(1)
external foo
c
c starting normalization is u = a x / ( c x + d ).
c follow this by something of the form v = u / ( copt u + dopt )
c (We optimize copt rather than c directly in a hope to deal
c with a better conditioned problem.)
c
dopt = 1 - copt*w(n)
lc = copt*a+dopt*c
ld = dopt*d
do 190 i=1,n-1
k = 2*i
ws = a*w(i)/(lc*w(i)+ld)
f(k-1) = real(ws-wt(i))
f(k) = aimag(ws-wt(i))
190 continue
end
subroutine sinter ( m, n, s, ijac, idj, jac, f, p )
c
c this routine defines the equation for the intersection
c of two circles, based on an expansion around the
c midpoint data.
c
real s(1), jac(idj,1), f(1), p(1)
integer n, ijac, idj
complex wm, wn, ws1, ws2, ci, e2
complex sexprl, e
real r1, r2, sr
ci= cmplx(0.e0,1.e0)
if(ijac.eq.1) goto 200
c
wm= cmplx(p(1),p(2))
wn= cmplx(p(3),p(4))
r1=p(5)
e=ci*s(1)/r1
e=sexprl(e)
e=s(1)*ci*e
ws1=wm+wn*e
wm= cmplx(p(6),p(7))
wn= cmplx(p(8),p(9))
r2=p(10)
e=ci*s(2)/r2
e=sexprl(e)
e=s(2)*ci*e
ws2=wm+wn*e
f(1)= real(ws1)- real(ws2)
f(2)=aimag(ws1)-aimag(ws2)
goto 900
c
200 continue
sr=s(1)/p(5)
e2= cmplx(p(3),p(4))* cmplx( cos(sr), sin(sr))
jac(1,1)=-aimag(e2)
jac(2,1)= real(e2)
sr=s(2)/p(10)
e2= cmplx(p(8),p(9))* cmplx( cos(sr), sin(sr))
jac(1,2)=aimag(e2)
jac(2,2)=- real(e2)
900 continue
return
end
subroutine smapc(c0, r, a, b, c, d, c3)
c
real r
complex c0, a, b, c, d, c3
c
c given the center c0 and radius r of a circle
c and the coefficients of a mobius transformation (a*z+b)/(c*z+d),
c compute the center c3 of the image of the circle.
c
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd, zk, ci, zero, wc(10), w0
c
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
c
real e, srinv
complex c1, c2, b1, b2
c
if (abs(c) .ge. eps) goto 10
c3 = (a*c0+b)/d
goto 40
10 c1 = c0+d/c
e = abs(c1)
if (e .ge. eps) goto 20
c3 = a/c
goto 40
20 if ( abs(e-r) .ge. 10*eps) goto 30
c3 = cmplx(srinv(0.0e0), 0.0e0)
goto 40
30 b1 = (1.0e0-r/e)*c1
b2 = (r/e+1.0e0)*c1
if(b1*b2*c.eq.0)write(6,*)"smapc",r,e,c1,b1,b2,c
c2 = .5e0*(1.e0/b1+1.e0/b2)
c3 = a/c+((b*c-a*d)/(c*c))*c2
40 continue
return
end
subroutine smtov(wm, wd)
c
complex wm(10), wd(2, 10)
c
c solve for vertices from midpoint information
c
common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun
integer maxfun
real radeps,odeeps,hybeps,epsfcn,dmax
c
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd, zk, ci, zero, wc(10), w0
c
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
c
common /sfail/ ifail
integer ifail
c
integer i, j, k, is, maxstp, iprint
real ws(2)
real tol, tr1(2), tr2(2,2), tr3(12), parm(10)
real u, v, marg(10), zarg(10)
real t, b1, b2, e
real theta, sarc, oarc
real det, a11, a12, a21, a22
complex w1, w2, ni, nj, zm, em, sexprl
external sinter
c
do 10 i = 1, n
zarg(i)=atan2(aimag(z(i)), real(z(i)))
if(zarg(i).lt.0) zarg(i) = zarg(i) + 2*pi
10 continue
marg(1) = zarg(1)/2.e0
is = n-1
do 20 i = 2, is
marg(i) = zarg(i-1)+(zarg(i)-zarg(i-1))/2.e0
20 continue
marg(n) = zarg(n-1)+(2.e0*pi-zarg(n-1))/2.e0
do 30 i = 1, n
zm = cmplx( cos(marg(i)), sin(marg(i)))
call scirc(zm, wm(i), wd(1, i), wd(2, i), wc(i), radc(i),
+ wdir(i))
c write(6,1025) i,wm(i),i,wd(1,i),i,wd(2,i),i,wc(i),
c + i,radc(i),i,wdir(i)
1025 format(
+ /' wm( ',i2,') = ',1p2e13.4
+ /' wd(1,',i2,') = ',1p2e13.4
+ /' wd(2,',i2,') = ',1p2e13.4
+ /' wc( ',i2,') = ',1p2e13.4
+ /' radc(',i2,') = ',1p1e13.4
+ /' wdir(',i2,') = ',i3 )
30 continue
c
c find intersection points
c
do 90 k = 1, n
i = k
j = mod(k, n)+1
if(abs(wc(i)-wc(j)).le.(radc(i)+radc(j))) goto 39
write(6,1039)
1039 format(' sides do not intersect!')
ifail=1
return
39 if (radc(j) .le. radc(i)) goto 40
i = j
j = k
40 if (radeps*radc(j) .le. 1.) goto 50
c
c two lines
c
b1 = ( real(wm(i)))* real(wd(1,i))
1 +(aimag(wm(i)))*aimag(wd(1,i))
b2 = ( real(wm(j)))* real(wd(1,j))
1 +(aimag(wm(j)))*aimag(wd(1,j))
a11 = real(wd(1, i))
a12 = aimag(wd(1, i))
a21 = real(wd(1, j))
a22 = aimag(wd(1, j))
det = (a11)*a22-(a12)*a21
w1=cmplx((b1*a22-a12*b2)/det,(-(a11*b1-b2*a21))/det)
w2=cmplx(1.e+15, 0.e0)
goto 70
50 if (radeps*radc(i) .le. 1.) goto 60
c
c circle and line
c
b1 = real(wd(1,i))
b2 = aimag(wd(1,i))
e = sqrt(b1**2+b2**2)
u = ( real(wd(1, i))*( real(wm(i)-wc(j)))
1 +aimag(wd(1, i))*(aimag(wm(i)-wc(j))))/e
v = sqrt((radc(j)-u)*(radc(j)+u))
w1 = wc(j)+ cmplx(u, v)*wd(1, i)/e
w2 = wc(j)+ cmplx(u, -v)*wd(1, i)/e
goto 70
c
c two circles
c
60 continue
if(radeps*radc(i).gt.radc(j) .or.
1 radeps*radc(j).gt.radc(i) ) then
write(6,1802)
1802 format(" circles with wildly differing radii")
ifail=3
return
endif
b1 = ( real(wc(j)))- real(wc(i))
b2 = (aimag(wc(j)))-aimag(wc(i))
e = sqrt(b1**2+b2**2)
b1 = radc(j)
b2 = radc(i)
u = (e**2-b2**2+b1**2)/(2.e0*e)
v = sqrt((b1-u)*(b1+u))
w1 = wc(j)+ cmplx(u, v)*(wc(i)-wc(j))/e
w2 = wc(j)+ cmplx(u, -v)*(wc(i)-wc(j))/e
c
c choose intersection point with correct interior angle
c
70 i = k
j = mod(k, n)+1
ni = w1-wc(i)
nj = w1-wc(j)
if(abs(ni).eq.0.0.or.abs(nj).eq.0.0) goto 801
ni = cmplx( float(wdir(i)), 0.e0)*ni/abs(ni)
nj = cmplx( float(wdir(j)), 0.e0)*nj/abs(nj)
theta=atan2(aimag(nj),real(nj))-atan2(aimag(ni),real(ni))
if ( abs(mod(alpha(k)-theta+9.e0*pi, 2.e0*pi)-pi)
1 .ge. 1.e0) go to 80
w(k) = w1
goto 85
80 w(k) = w2
85 continue
zm= cmplx( cos(marg(i)), sin(marg(i)))
ni=wdir(i)*wd(1,i)*zm
ni=ni/abs(ni)
zm= cmplx( cos(marg(j)), sin(marg(j)))
nj=wdir(j)*wd(1,j)*zm
nj=nj/abs(nj)
c write(6,1031)i,ni,j,nj,w(i)
1031 format(
+ /' snewt2'
+ /' i=',i2
+ ,' ni=',0p2f7.3
+ /' j=',i2
+ ,' nj=',0p2f7.3
+ /' w before =',0p2e12.4 )
parm(1)= real(wm(i))
parm(2)=aimag(wm(i))
parm(3)= real(ni)
parm(4)=aimag(ni)
parm(5)=radc(i)
parm(6)= real(wm(j))
parm(7)=aimag(wm(j))
parm(8)= real(nj)
parm(9)=aimag(nj)
parm(10)=radc(j)
tol=5*eps*abs(w(k))
maxstp=10
iprint=0
c**** 20 Apr 85 added wdir
c write(6,1803) w(k),wm(i),wc(i),wm(j),wc(j)
c1803 format('for',1p10e10.2)
c ws(1)=sarc(w(k),wm(i),wc(i),1)
c ws(2)=sarc(w(k),wm(j),wc(j),1)
c write(6,1804) ws(1), ws(2)
c1804 format('new',1p10e10.2)
ws(1)=oarc(w(k),wm(i),wc(i))
ws(2)=oarc(w(k),wm(j),wc(j))
c write(6,1805) ws(1), ws(2)
c1805 format('old',1p10e10.2)
call snewt2(2,2,ws,tr1,tol,maxstp,1.0e0,iprint,
+ 2,tr2,sinter,parm,tr3)
em=ci*ws(1)*sexprl(ci*ws(1)/radc(i))
w(k)=wm(i)+ni*em
c write(6,1035)w(k)
1035 format(' w after =',2e12.4 )
90 continue
return
801 continue
write(6,1801)
1801 format(' smtov: entire side mapped onto a point')
ifail=2
return
end
subroutine snewt2(m,n,x,f,tol,maxstp,dmax,
+ iprint,idj,jac,fun,parm,w)
c
integer m,n,idj,maxstp,iprint
real x(1),f(1),jac(idj,1)
real dmax
real w(1)
real parm(1)
external fun
c
c This is a simple Newton code for solving a system of nonlinear
c equations. If the number of equations m exceeds the number of
c unknowns n then a nonlinear least squares solution is attempted.
c
c dmax is the maximum 2-norm of any step that will be attempted.
c
c The user must supply a subroutine fun, declared external in the
c calling program , with the following callingsequence:
c
c fun(m,n,x,ijac,idj,jac,f,parm)
c
c if ijac=0 : f should return the function value at x.
c if ijac=1 : jac should return the jacobian Df(i)/Dx(j).
c In this case f is input and already contains
c the function value at x. (from a previous call)
c
c This code will return when the 2-norm of f is less than tol or
c when maxstp steps have been tried. (Each involving one jacobian
c evaluation and at most 3 function evaluations.)
c
c The code may also give up and in this case some information
c about the function is written out.
c
c w is a workspace containing at least 5*n+m elements.
c parm is a real array in which the user can pass information
c to fun. It is never referenced by Newton.
c
c Local variables.
c
integer i1,i2,i3,i4,i5
real fnrm
c
if(maxstp.lt.1) go to 100
if(m.lt.n) go to 200
c
c assign workspace.
c
i1 = n+1
i2 = i1+m
i3 = i2+n
i4 = i3+n
i5 = i4+n
c
nfun = 0
njac = 0
c
call snwt2(m,n,x,f,tol,maxstp,dmax,iprint,
+ idj,jac,fun,parm,nfun,njac,fnrm,
+ w(1),w(i1),w(i2),w(i3),w(i4),w(i5))
c
if(iprint.le.0) return
write(6,1) fnrm
write(6,2) nfun
write(6,3) njac
return
100 write(6,4)
return
200 write(6,5)
return
1 format(1x,'on exit from newton the norm of f is,'e12.3)
2 format(1x,'number of function evaluations were ,'i6)
3 format(1x,'number of jacobian evaluations were ,'i6)
4 format(1x,'newton returned because maxstp was ,'i6)
5 format(1x,'newton returned because m is less than n, m,n=',2i4)
end
subroutine snwt2(m,n,x,f,tol,maxstp,dmax,iprint,
+ idj,jac,fun,parm,nfun,njac,fnrm,
+ xnew,fnew,stepn,qraux,jpvt,stepg)
c
integer m,n,maxstp,iprint,idj,nfun,njac
real tol,dmax,fnrm
integer jpvt(1)
real x(1),f(1)
real jac(idj,1)
real xnew(1),fnew(1),stepn(1),qraux(1),stepg(1)
real parm(1)
c
c this routine must always be called from the
c driver routine newtn.
c
c local variables.
c
integer i,j,it,info
real fnrm1,fnrm2,fnrm3,fmin
real dec1,dec2
real a,b,c,t
real alfa
real sdot,snrm2
real tmp(10,10),xtmp(10),ytmp(10),ztmp(10)
c
c dec1 is used to test for decrease.
c dec2 is used to test for decrease if dec1 cannot be satisfied.
c alfa is steplenght parameter.
c
dec1 = .2e0
dec2 = .9e0
call scopy(n,x,1,xnew,1)
nfun = nfun + 1
call fun(m,n,x,0,idj,jac,f,parm)
fnrm = snrm2(m,f,1)
if(fnrm.lt.tol) go to 600
fmin = fnrm
c
c start main iteration loop.
c
do 400 it=1,maxstp
if(iprint.le.0) go to 210
write(6,10) it,fnrm
write(6,11)(x(i),i=1,n)
write(6,12)(f(i),i=1,m)
210 njac = njac + 1
call fun(m,n,x,1,idj,jac,f,parm)
c
c if iprint is 3 or larger then diagnostic info
c about the jacobian can be (computed and) printed
c here.
c
if(iprint.le.2) go to 220
write(6,18)
write(6,19)((jac(i,j),j=1,n),i=1,m)
do 215 i=1,n
do 215 j=1,n
tmp(i,j)=jac(i,j)
215 continue
do 216 i=1,n
ytmp(i)=f(i)
216 continue
call sgeco(tmp,10,n,xtmp,rcond,ztmp)
call sgesl(tmp,10,n,xtmp,ytmp,0)
write(6,719) rcond
write(6,718)(ztmp(i),i=1,n)
write(6,718)(ytmp(i),i=1,n)
719 format(1x,'rcond=',e12.3)
718 format(1x,4e12.3)
220 do 230 i = 1,n
jpvt(i) = 0
230 continue
call sqrdc(jac,idj,m,n,qraux,jpvt,stepg,1)
call sqrsl(jac,idj,m,n,qraux,f,f,f,stepg,f,f,100,info)
if(info.ne.0) go to 520
do 240 i = 1,n
stepn(jpvt(i)) = -stepg(i)
240 continue
if(iprint.le.2) go to 250
write(6,20)(stepn(i),i=1,n)
c
c test if newton step is too large.
c
250 fnrm1 = snrm2(n,stepn,1)
if(fnrm1.lt.dmax) go to 258
if(iprint.ge.3) write(6,22) fnrm1,dmax
c
c the gradient of the squared norm of f is computed and
c normalized to lenght dmax. The step is taken to be
c between the newton step and this
c gradient direction. (See below for details.)
c
do 256 j = 1,n
t = 0.e0
do 255 i = 1,m
t = t - f(i)*jac(i,j)
255 continue
stepg(j) = t
256 continue
c t = dmax/fnrm1
call sscal(n,1.e0/fnrm1,stepn,1)
fnrm1 = snrm2(n,stepg,1)
call sscal(n,-1.e0/fnrm1,stepg,1)
fnrm1 = sdot(n,stepn,1,stepg,1)
if(iprint.ge.2) write(6,17) fnrm1
c do 257 i = 1,n
c stepn(i) = t*stepn(i) + (1.e0-t)*stepg(i)
c 257 continue
call sscal(n,dmax,stepn,1)
c
c try a unit step.
c
258 call saxpy(n,1.e0,stepn,1,xnew,1)
if(iprint.ge.3) write(6,21)(xnew(i),i=1,n)
nfun = nfun + 1
call fun(m,n,xnew,0,idj,jac,fnew,parm)
fnrm1 = snrm2(m,fnew,1)
if(iprint.le.1) go to 260
write(6,13) fnrm1
write(6,11) (xnew(i),i=1,n)
write(6,12) (fnew(i),i=1,m)
260 fmin = min(fmin,fnrm1)
if(fnrm1.lt.dec1*fnrm) go to 350
if(fnrm1.gt.fmin) go to 270
alfa = 1.e0
call scopy(m,fnew,1,f,1)
c
c try half a unit step.
c
270 call saxpy(n,-.5e0,stepn,1,xnew,1)
nfun = nfun + 1
call fun(m,n,xnew,0,idj,jac,fnew,parm)
fnrm2 = snrm2(m,fnew,1)
if(iprint.le.1) go to 280
write(6,14) fnrm2
write(6,11)(xnew(i),i=1,n)
write(6,12)(fnew(i),i=1,m)
280 fmin = min(fmin,fnrm2)
if(fnrm2.lt.dec1*fnrm) go to 350
if(fnrm2.gt.fmin) go to 290
alfa = .5e0
call scopy(m,fnew,1,f,1)
c
c compute the quadratic passing through fnrm, fnrm2 and fnrm1.
c
290 a = 2.e0*(fnrm+fnrm1-2.e0*fnrm2)
b = 4.e0*fnrm2-3.e0*fnrm-fnrm1
c = fnrm
if(a.le.0.e0) go to 310
t = -b/(2.e0*a)-.5e0
c
c try a step to the minimum of a quadratic model.
c
call saxpy(n,t,stepn,1,xnew,1)
t = t + .5e0
nfun = nfun + 1
call fun(m,n,xnew,0,idj,jac,fnew,parm)
fnrm3 = snrm2(m,fnew,1)
if(iprint.le.1) go to 300
write(6,15) a,b,c
write(6,16) fnrm3,t
write(6,11)(xnew(i),i=1,n)
write(6,12)(fnew(i),i=1,m)
300 fmin = min(fmin,fnrm3)
if(fnrm3.lt..5e0*(dec1+dec2)*fnrm) go to 350
if(fnrm3.gt.fmin) go to 310
alfa = t
call scopy(m,fnew,1,f,1)
310 if(fmin.gt.dec2*fnrm) go to 530
call scopy(n,x,1,xnew,1)
call saxpy(n,alfa,stepn,1,xnew,1)
go to 360
c
c update and prepare for a new iteration.
c
350 call scopy(m,fnew,1,f,1)
360 call scopy(n,xnew,1,x,1)
fnrm = fmin
if(fnrm.lt.tol) go to 600
400 continue
write(6,1)
return
520 write(6,3)
return
530 write(6,4) fnrm
write(6,11)(x(i),i=1,n)
write(6,7)(stepn(i),i=1,n)
write(6,5) fnrm,fnrm2,fnrm1
write(6,6) alfa,fmin
write(6,11)(xnew(i),i=1,n)
write(6,12)(f(i),i=1,m)
return
600 continue
return
1 format(1x,'maxstp iterations in newtn')
3 format(1x,'jacobian has rank less than n')
4 format(1x,'Newton gave up, norm of f is:'e12.4)
5 format(1x,'Norm of f in x+alfa*stepn , alfa=0,1/2,1, is :'3e12.3)
6 format(1x,'With step= ',e12.4,1x,',the norm is :',e12.3)
7 format(1x,'s= ',1x,5e14.6,4(/5x,5e14.6))
10 format(/1x,'At the start of iteration',i3,4x,'Norm of f= ',e14.6)
11 format(1x,'x= ',1x,5e14.6,4(/5x,5e14.6))
12 format(1x,'f= ',1x,5e14.6,4(/5x,5e14.6))
13 format(1x,'Unit step,',1x,'Norm of f=',e12.4)
14 format(1x,'Half step,',1x,'Norm of f=',e12.4)
15 format(1x,'Quadratic model coeff a,b,c=',1x,3e12.4)
16 format(1x,'Min step,',1x,'Norm of f=',e12.4,2x,'alfa=',e12.4)
17 format(1x,'angle between newton and gradient direction,',e12.3)
18 format(1x,'jacobian matrix.')
19 format(10(5x,5e12.3/))
20 format(1x,'newton step,'/,5(5x,5e12.3/))
21 format(1x,'newton step gives x=,'/,5(5x,5e12.3/))
22 format(1x,'step too large, norm of step,dmax= ',2e12.3)
end
subroutine sray(theta,y0,y1,y2,ntrace,trace)
c
c trace radial ray in unit disk to unnormalized polygon
c
integer ntrace
real theta
complex y0, y1, y2, trace(ntrace)
c
common /initc/ a, b, c, d, winit1, winit2
complex a, b, c, d, winit1, winit2
c
common /epscom/ radeps,odeeps,hybeps,epsfcn,dmax,maxfun
integer maxfun
real radeps,odeeps,hybeps,epsfcn,dmax
c
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd, zk, ci, zero, wc(10), w0
c
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
c
common /savex/ xxx
real xxx(10)
common/node/ nfun
integer nfun
c
external sfode, sfode4
integer iflag, iwork(50), odeit, odemit, step
real y(6), relerr, abserr, r1, r2, rwork(600)
real yeps(6), ymag(6)
common /sfail/ ifail
integer ifail
integer istate, jt, index, order, limits(7)
real dt, dtmax
data odemit /10/
c
zd = cmplx( cos(theta), sin(theta))
y(1) = 0.e0
y(2) = 0.e0
y(3) = real(winit1*zd)
y(4) = aimag(winit1*zd)
y(5) = real(winit2*zd**2)
y(6) = aimag(winit2*zd**2)
r1 = 0.e0
relerr = odeeps
abserr = 1.e-5*relerr
iflag = 1
odeit = 1
trace(1) = cmplx(y(1),y(2))
do 50 step=2,ntrace
r2 = (step-1.0e0)/(ntrace-1.0e0)
c
c istate=1
c jt=2
c call lsoda(sfode4,6,y,r1,r2,1,relerr,abserr,1,istate,
c $ 0,rwork,200,iwork,26,sfode4,jt)
c if(istate.ne.2) write(6,*)"lsoda returned",istate
c...
c index=0
c dt=.001
c dtmax=.1
c do 10 i=1,6
c yeps(i) = relerr
c 10 continue
c 11 call stride(6,sfode,sfode,r1,r2,y,ymag,index,dt,dtmax,
c $ order,yeps,101,limits,iwork,rwork)
c if(index.ge.1) goto 11
cc...
c goto 40
21 call ode(sfode,6,y,r1,r2,relerr,abserr,iflag,rwork,iwork)
if(iflag.eq.2) goto 40
if(iflag.ne.3) goto 22
write(6,1021) relerr,abserr
1021 format(" ode increased relerr,abserr to ",1p2e11.2)
odeit=odeit+1
if(odeit.le.odemit)goto 21
ifail=1
return
c
22 if(iflag.ne.4) goto 23
write(6,1022)
1022 format(" ode complained that it had used more than 500 steps")
odeit=odeit+1
if(odeit.le.odemit)goto 21
ifail=1
return
c
23 if(iflag.ne.5) goto 24
write(6,1023) r1, r2
1023 format(" ode complained equations are stiff r=",2f8.4)
c------
m = 2*n-4
write(6,8024) (i,xxx(i),i=1,m)
8024 format(' x(',i2,')=',f22.16)
y(1) = 0.e0
y(2) = 0.e0
y(3) = real(winit1*zd)
y(4) = aimag(winit1*zd)
y(5) = real(winit2*zd**2)
y(6) = aimag(winit2*zd**2)
do 8023 i=1,6
write(6,*)" y(",i,")=",y(i)
8023 continue
call sfode(-1,y,yp)
stop
c------
odeit=odeit+1
if(odeit.le.odemit)goto 21
ifail=1
return
c
24 write(6,1) iflag
1 format(1x,'ode returned iflag = ',1x,i8/)
ifail=1
return
c
40 continue
trace(step) = cmplx(y(1),y(2))
50 continue
c
y0 = cmplx(y(1), y(2))
y1 = cmplx(y(3), y(4))/zd
y2 = cmplx(y(5), y(6))/zd**2
return
end
function srinv(x)
c
real x
c
if ( abs(x) .ge. 1.e-15) goto 1
srinv = sign(1.e+15, x)
goto 2
1 srinv = 1.0e0/x
2 continue
return
end
subroutine sstpar(x)
c
real x(1)
c
c given guesses of vertex preimages and n-3 of the
c auxiliary parameters, solve for remaining 3 parameters
c
common /capcm1/ wt, cnr, wtc, z, cc, w, wc, w0, zd, zk, ci, zero
complex wt(10), cnr, wtc(10), z(10), cc(10), w(10)
complex zd, zk, ci, zero, wc(10), w0
c
common /capcm2/ alpha, gamma, beta, s, pi, eps, n, it, wtdir,
1 wdir, iv, radt, radc, arct, arcc
integer n, it, wtdir(10), wdir(10), iv(1)
real s(6)
real alpha(10), gamma(10), beta(10), pi, eps
real radt(10), radc(10), arct(10), arcc(10)
c
common /savex/ xxx
real xxx(10)
integer i, m
real a(10)
call strans(x, z, n)
m = 2*n-4
do 20 i = 1, n
cc(i)= ci*(zk+z(i))/(zk-z(i))
a(i) = real(cc(i))
xxx(i) = x(i)
20 continue
do 30 i = n, m
beta(i-n+1) = x(i)
30 continue
call selim(n, a, beta, gamma)
return
end
subroutine strans(x,z,n)
integer n
real x(1), a(25)
complex z(n)
if(n.gt.25)then
write(6,*)'strans a dimension too samll'
stop
end if
call stransa(x,z,a,n)
return
end
subroutine stransa(x,z,a,n)
c
c given optimization parameters x, back transform to
c get angles a and points z on unit circle
c
c x(i)=log((teta(i)-teta(i-1))/(teta(i+1)-teta(i)))
c for i=1,2.. n-1.
c teta(0)=0,teta(n)=twopi.
c z(i)=exp(%i*teta(i)).
c
integer n
real x(1), a(n)
complex z(n)
integer i,nm
real y,sum,twopi
twopi=8.e0*atan(1.e0)
y=1.e0
sum=y
nm=n-1
do 10 i=1,nm
y=y/exp(x(i))
sum=sum + y
10 continue
y=twopi/sum
sum=y
a(1)=sum
z(1)= cmplx( cos(y), sin(y))
do 20 i=2,nm
y=y/exp(x(i-1))
sum=sum+y
a(i)=sum
z(i)= cmplx( cos(sum), sin(sum))
20 continue
a(n)=0.e0
z(n)= cmplx(1.e0,0.e0)
return
end
c From jcb Thu Feb 6 15:16 EST 1986
SUBROUTINE STRIDE(N,F,JAC,X,XPRINT,Y,YMAG,INDEX,STP,STPMAX,M,EPS,
+ METHOD,LIMITS,INTGRS,REALS)
COMMENT OOO:
C
C S T R I D E
C
C STABLE RUNGE-KUTTA INTEGRATOR FOR DIFFERENTIAL EQUATIONS
C
C J C BUTCHER, K BURRAGE, F H CHIPMAN
C
C MARCH 1979
C
C
C THIS VERSION OF STRIDE CONTAINS THE CHANGES THAT WERE SPECIFIED
C IN F.CHIPMANS LETTER TO P.W.GAFFNEY ON SEPTEMBER 21,1979.
C
COMMENT 020: DECLARATIONS
C
DIMENSION Y(1),YMAG(1),EPS(1),LIMITS(1),INTGRS(1),REALS(1)
INTEGER YTEMP,YSAVE,FTOTAL,FVECT,YVECT,FSPACE,YSPACE,EFACTS,
+T,TINV,XI,DERIV,YRES,W,YLAST,FACTS
LOGICAL CONVGD,DIVGD,EXCESS,REVAL,OK
C
COMMENT 040: IF NOT FIRST CALL UNPACK REALS AND INTGRS
C
IF(INDEX.EQ.0)GO TO 120
GO TO 080
C
COMMENT 060: STATEMENTS REPRESENTING THE 'PAUSE' ROUTINE
C
060 IF(INDEX .GT. 1) LIMITS(INDEX)=LIMITS(INDEX)-1
IF (INDEX .GT. 2) LIMITS(1)=LIMITS(1)-1
IF (INDEX.LT.3) GO TO 070
IF (INTGRS(LIMMEM+INDEX).GT.0 .AND.
+ LIMITS(INDEX)+INTGRS(LIMMEM+INDEX).LT.0) GOTO 070
GOTO 100
C THEN
070 IF(INDEX.GE.3)LIMITS(INDEX)=INTGRS(LIMMEM+INDEX)
IF(INDEX.EQ.1)RETURN
INTGRS(1) = MLAST
INTGRS(3) = MSTATE
INTGRS(4) = MTRAN
INTGRS(5) = NJAC
INTGRS(6) = NSTEP
INTGRS(7) = 0
IF(REVAL)INTGRS(7) = 1
REALS(2) = EGROW
REALS(3) = ERREST
REALS(4) = H
REALS(5) = HFACTS
REALS(6) = HLAST
REALS(7) = XRES
REALS(8) = XLAST
REALS(9) = XJAC
RETURN
080 IF(INDEX.GT.1)GO TO 085
MAXM = M
IF(MAXM.GT.15)MAXM=15
INTGRS(2) = MAXM
085 URND = REALS(1)
MAXM = INTGRS(2)
MAXI=7
LIMMEM = MAXI + MAXM + 1
IPERM = LIMMEM + 7
EFACTS = 9
FVECT = EFACTS + INTGRS(MAXI+MAXM+1)
YVECT = FVECT + MAXM
YSPACE = YVECT + MAXM - N
FSPACE = YSPACE + (MAXM + 2)*N
T = FSPACE + (MAXM + 2)*N + N - MAXM
TINV = T + MAXM**2
XI = TINV + MAXM**2 + MAXM
DERIV = XI + MAXM*(MAXM + 1)/2
YRES = DERIV + N
YTEMP = YRES + N
YSAVE = YTEMP + N
W = YSAVE + N
FTOTAL = W + N
YLAST = FTOTAL + N
FACTS = YLAST + N - N
MERPST=METHOD / 100
MTYPE=METHOD / 10
MAXITS=METHOD-10*MTYPE
MTYPE=MTYPE-10*MERPST
IF(INDEX.NE.1)GO TO 090
GO TO 095
C THEN
090 MLAST = INTGRS(1)
MSTATE = INTGRS(3)
MTRAN = INTGRS(4)
NJAC = INTGRS(5)
NSTEP = INTGRS(6)
REVAL = INTGRS(7).EQ.1
EGROW = REALS(2)
ERREST = REALS(3)
H = REALS(4)
HFACTS = REALS(5)
HLAST = REALS(6)
XRES = REALS(7)
XLAST = REALS(8)
XJAC = REALS(9)
C FI
095 CONTINUE
C FI
100 DO 105 I=1,7
IF(LIMITS(I).GE.0)GO TO 101
GO TO 105
C THEN
101 INTGRS(LIMMEM+I)=LIMITS(I)
LIMITS(I) = -1
C FI
105 CONTINUE
C OD
DO 110 I=1,2
110 IF(INTGRS(LIMMEM+I).GT.0.AND.LIMITS(I)+
+ INTGRS(LIMMEM+I).LT.0)INDEX=0
C OD
IF ( INDEX .GT. 3 ) STP=0
IF(INDEX.EQ.0)RETURN
GO TO (140,860,915,760,760,820,885),INDEX
C
COMMENT 120: INITIALIZE CONSTANTS AND STARTING VALUES
C
120 MAXI=7
TEMP = 1.E0
URND=1.E0
125 URND = URND/2.E0
IF (1.E0 + URND/2.E0 .GT. TEMP) GO TO 125
REALS(1) = URND
INTGRS(MAXI + 1) = 0
INTGRS(MAXI + 2) = INTGRS(MAXI + 1) + 1
INTGRS(MAXI + 3) = INTGRS(MAXI + 2) + 2
INTGRS(MAXI + 4) = INTGRS(MAXI + 3) + 2
INTGRS(MAXI + 5) = INTGRS(MAXI + 4) + 3
INTGRS(MAXI + 6) = INTGRS(MAXI + 5) + 3
INTGRS(MAXI + 7) = INTGRS(MAXI + 6) + 4
INTGRS(MAXI + 8) = INTGRS(MAXI + 7) + 4
INTGRS(MAXI + 9) = INTGRS(MAXI + 8) + 4
INTGRS(MAXI + 10) = INTGRS(MAXI + 9) + 5
INTGRS(MAXI + 11) = INTGRS(MAXI + 10) + 5
INTGRS(MAXI + 12) = INTGRS(MAXI + 11) + 5
INTGRS(MAXI + 13) = INTGRS(MAXI + 12) + 6
INTGRS(MAXI + 14) = INTGRS(MAXI + 13) + 6
INTGRS(MAXI + 15) = INTGRS(MAXI + 14) + 6
INTGRS(MAXI + 16) = INTGRS(MAXI + 15) + 6
DO 130 I=1 , N
130 YMAG(I)=0.E0
C OD
M = 10
DO 135 I=1,7
135 LIMITS(I) = 0
LIMITS(2) = 1
STPMAX=0.E0
STP=0.E0
INDEX=1
GO TO 060
C GO TO AND RETURN FROM PAUSE ROUTINE
140 H=STP
XLAST=X
DO 145 I=1 , N
145 IF( ABS(Y(I)).GT.YMAG(I) ) YMAG(I)=ABS(Y(I))
C OD
XRES=0.E0
DO 150 J=1 , N
150 REALS(YRES+J)=0.E0
C OD
IF(H.EQ.0.E0) GOTO 155
GOTO 175
C THEN
155 ENORM=0.E0
DO 160 I=1 , N
160 IF( YMAG(I).GT.EPS(I)*ENORM ) ENORM=YMAG(I)/EPS(I)
C OD
IF(ENORM.NE.0.E0) GOTO 165
GOTO 175
C THEN
165 CALL F(X,Y,REALS(YTEMP+1))
DO 170 I = 1 , N
170 IF(REALS(YTEMP+I).NE.0.E0.AND. (H.EQ.0.E0
+ .OR.H*ABS(REALS(YTEMP+I)).GT.1.25E0*EPS(I)*ENORM-YMAG(I)))
+ H=(1.25E0*EPS(I)*ENORM-YMAG(I))/ABS(REALS(YTEMP+I))
C OD
C FI
C FI
175 IF( H.EQ.0.E0 ) H=XPRINT-X
IF(ABS(H).GT.STPMAX .AND. STPMAX.GT.0.E0) H=STPMAX
IF(H*(XPRINT-X).LT.0.E0) H=-H
MLAST = 0
M=1
REVAL = .TRUE.
NJAC=0
STP=0.E0
MSTATE = 0
MTRAN=0
HFACTS=0.E0
HLAST=0.E0
ERREST=1.E0
XJAC=X
C
COMMENT 180: COMPUTE LAGUERRE ZEROS AND RELATED CONSTANTS
C
REALS(XI+1) =1
REALS(EFACTS+1)=0.5E0
DO 205 M=2 , MAXM
DO 205 I=1 , M
K= M * (M-1) / 2
IF (I.EQ.1) GO TO 185
IF (I.EQ.M) GO TO 190
TEMP= (REALS(XI+K-M+I) + REALS(XI+K-M+I+1))/2.E0
GO TO 195
185 TEMP= REALS(XI+K-M+2) /2.E0
GO TO 195
190 TEMP= REALS(XI+K)*(1.E0+1.E0/(M-1.9E0))
195 P= 1
PPRIME= 0.E0
DO 200 J= 1 , M
PPRIME= PPRIME - P
200 P= P + TEMP * PPRIME/J
C OD
DELTA=P/PPRIME
TEMP=TEMP-DELTA
IF( (DELTA/TEMP) ** 2 .GE. .01E0 * URND )GO TO 195
REALS(XI+K+I) = TEMP
205 IF( I.LE.INTGRS(MAXI+M+1)-INTGRS(MAXI+M) )
+ REALS(EFACTS+INTGRS(MAXI+M)+I)=ABS(PPRIME*TEMP)/(M+1)
C OD
C OD
C
COMMENT 220: PERFORM STEPS
C
M = 1
NSTEP=0
C DO
220 NSTEP=NSTEP+1
IF( H*(XPRINT-X).LT.0.E0 ) GO TO 225
GO TO 230
C THEN
225 MSTATE=0
H=-H
MLAST=0
C FI
230 MROOTS=M*(M-1) / 2
REVAL = REVAL.AND.MTYPE .GT. 1
TEMP=HFACTS/H
IF(.NOT.REVAL.AND.MTYPE.GT.1)
+ REVAL=TEMP.GT.1.5625E0.OR.TEMP.LT.0.64E0.OR. NSTEP.GE.NJAC + 20
IF(MTYPE.LE.1)HFACTS=0.E0
IF (REVAL) GOTO 240
GOTO 360
C THEN
C
COMMENT 240: RE-EVALUATE JACOBIAN IF APPROPRIATE
C
240 IF( MTYPE.EQ.2 ) GO TO 241
GO TO 245
C THEN
241 CALL JAC(X,Y,REALS(FACTS+N+1))
GO TO 265
C ELSE
245 SQURND=SQRT(URND)
CALL F(X,Y,REALS(W+1))
D = 0.E0
DO 250 I = 1 , N
250 D = D + ABS(REALS(W+I))/EPS(I)
C OD
D = D * ABS(H) * URND * 1.E-3
DO 260 J = 1 , N
DELTA = D * EPS(J) + SQURND * YMAG(J)
TEMP = Y(J)
Y(J) = TEMP + DELTA
DELTA = Y(J) - TEMP
CALL F(X,Y,REALS(DERIV+1))
DO 255 I = 1 , N
REALS(FACTS+N*J+I) = 0.E0
255 IF(DELTA .NE. 0.E0) REALS(FACTS+N*J+I)=
+ (REALS(DERIV + I)-REALS(W + I))/DELTA
260 Y(J) = TEMP
C OD
C OD
C FI
265 XJAC=X
NJAC=NSTEP
REVAL=.FALSE.
C
COMMENT 280: FACTORIZE ( IDENTITY - H * JACOBIAN )
C
HFACTS= H
DO 290 I= 1 , N
DO 285 J= 1 , N
285 REALS(FACTS+N*J+I)=-H*REALS(FACTS+N*J+I)
C OD
REALS(FACTS+N*I+I)=1+REALS(FACTS+N*I+I)
290 INTGRS(IPERM+I)= I
C OD
DO 350 K= 1 , N
IF(.NOT.REVAL) GOTO 291
GOTO 350
C THEN
291 PIVOT= 0.E0
IPIVOT= 0
IF(K .GT. 1) GOTO 295
GOTO 330
C THEN
295 DO 325 I =1,N
TEMP= REALS(FACTS+N*K+I)
IF(I .LT. K) GOTO 296
GOTO 310
C THEN
296 IF(I.EQ.1)GO TO 305
IM1 = I-1
DO 300 J=1 , IM1
300 TEMP= TEMP-REALS(FACTS+N*J+I)*REALS(FACTS+N*K+J)
C OD
305 REALS(FACTS+N*K+I)= TEMP/REALS(FACTS+N*I+I)
GO TO 325
C ELSE
310 KM1 = K-1
IF(K.EQ.1)GO TO 320
DO 315 J=1 , KM1
315 TEMP= TEMP-REALS(FACTS+N*J+I)*REALS(FACTS+N*K+J)
C OD
320 REALS(FACTS+N*K+I)= TEMP
IF( ABS(TEMP/EPS(INTGRS(IPERM+I))) .GT. PIVOT )
+ GOTO 321
GOTO 325
C THEN
321 PIVOT= ABS(TEMP/EPS(INTGRS(IPERM+I)))
IPIVOT= I
C FI
C FI
325 CONTINUE
C OD
GO TO 350
C ELSE
330 DO 335 I= 1 , N
IF( ABS(REALS(FACTS+N+I)/EPS(I)) .GT. PIVOT )
+ GOTO 331
GOTO 335
C THEN
331 PIVOT= ABS(REALS(FACTS+N+I)/EPS(I))
IPIVOT= I
335 CONTINUE
C FI
C OD
IF( IPIVOT .EQ. 0 ) GOTO 336
GOTO 340
C THEN
336 H = 0.9E0*H
MSTATE = 1
GO TO 240
C ELSE
340 IF( IPIVOT .NE. K )GO TO 341
GOTO 350
C THEN
341 J= INTGRS(IPERM+IPIVOT)
INTGRS(IPERM+IPIVOT)= INTGRS(IPERM+K)
INTGRS(IPERM+K)= J
DO 345 J= 1 , N
TEMP = REALS(FACTS+N*J+IPIVOT)
REALS(FACTS+N*J+IPIVOT) = REALS(FACTS+N*J+K)
345 REALS(FACTS+N*J+K) = TEMP
C OD
C FI
C FI
C FI
350 CONTINUE
C FI
C OD
C FI
C
COMMENT 360: COMPUTE STARTING VALUES FOR ITERATIONS
C
360 IF( STP .NE. 0.E0 .OR. H .NE. HLAST )GO TO 361
GOTO 440
C THEN
361 IF( MLAST.GT.M ) MLAST=M
IF( MLAST.EQ.0 )GO TO 365
GOTO 375
C THEN
365 MP1 = M+1
DO 370 I=1 , MP1
DO 370 J=1 , N
370 REALS(YSPACE+N*I+J)=0.E0
C OD
C OD
GO TO 440
C ELSE
375 MP1 = M+1
DO 425 KP1 = 1,MP1
K = KP1 - 1
D = 0.E0
IF(K.NE.0)D=H*REALS(XI+MROOTS+K)
D = (STP + D)/HLAST
P=1.E0
PPRIME=0.E0
TEMP=0.E0
IF( K.EQ.0 )GO TO 380
GOTO 390
C THEN
380 DO 385 J=1,N
REALS(YSPACE+N*M+N+J)=0.E0
385 REALS(YTEMP+J)=0.E0
C OD
GO TO 400
C ELSE
390 DO 395 J=1 , N
395 REALS(YSPACE+N*K+J)=-REALS(YTEMP+J)
C OD
C FI
400 DO 425 I=1 , MLAST
PPRIME=PPRIME-P
TEMP=(D/I)*PPRIME
P=P+TEMP
IF( K.EQ.0 )GO TO 405
GOTO 415
C THEN
405 DO 410 J = 1,N
410 REALS(YTEMP+J)=REALS(YTEMP+J)-TEMP *
+ REALS(FSPACE+N*I+J)
C OD
GO TO 425
C ELSE
415 DO 420 J=1 , N
420 REALS(YSPACE+N*K+J)=
+ REALS(YSPACE+N*K+J)-TEMP*REALS(FSPACE+N*I+J)
C OD
C FI
C OD
C OD
C FI
C FI
425 CONTINUE
C
COMMENT 440: GENERATE TRANSFORMATION MATRICES T AND TINV
C WHEN NECESSARY
C
440 IF( M .NE. MTRAN )GO TO 441
GOTO 460
C THEN
441 MTRAN = M
DO 455 I= 1 , M
P= 1.E0
PPRIME= 0.E0
TEMP= REALS(XI+MROOTS + I)
REALS(T+MAXM+I)= P
IF(M.EQ.1)GO TO 450
DO 445 J= 2 , M
PPRIME= PPRIME - P
P= P + TEMP*PPRIME/(J-1)
445 REALS(T+MAXM*J+I)= P
C OD
450 DO 455 K= 1 , M
455 REALS(TINV+MAXM*I+K)= TEMP*REALS(T+MAXM*K+I)/(M*P)**2
C OD
C OD
C FI
C
COMMENT 460: CARRY OUT ITERATIONS
C
460 NITS = 0
CONVGD = .FALSE.
DIVGD = .FALSE.
EXCESS = .FALSE.
465 NITS = NITS + 1
IF(.NOT.(CONVGD.OR.DIVGD.OR.EXCESS))GO TO 480
GOTO 740
C THEN
C
COMMENT 480: EVALUATE DERIVATIVES AND TRANSFORM USING TINV
C
480 ISAVE=INTGRS(MAXI+M+1)-INTGRS(MAXI+M)
DO 485 J =1 , N
485 REALS(YSAVE+J)=REALS(YSPACE+N*ISAVE+J)
C OD
MP1 = M+1
DO 495 I=1 , MP1
DO 490 J=1 , N
490 REALS(YTEMP+J)= Y(J) + REALS(YSPACE+N*I+J)
C OD
TEMP = X
IF (I .LE. M) TEMP = X+H*REALS(XI+MROOTS+I)
CALL F(TEMP,REALS(YTEMP+1),REALS(DERIV+1))
DO 495 K=1 , N
495 REALS(FSPACE+N*I+K)= H*REALS(DERIV+K)
C OD
C OD
DO 510 J=1 , N
DO 505 I=1 , M
TEMP=0.E0
TEMPF=0.E0
DO 500 K=1 , M
TEMP=TEMP + REALS(TINV+MAXM*K+I)*REALS(YSPACE+N*K+J)
500 TEMPF=TEMPF+ REALS(TINV+MAXM*K+I)*REALS(FSPACE+N*K+J)
C OD
REALS(YVECT+I)= TEMP
505 REALS(FVECT+I)= TEMPF
C OD
DO 510 I=1 , M
REALS(YSPACE+N*I+J)= REALS(YVECT+I)
510 REALS(FSPACE+N*I+J)= REALS(FVECT+I)
C OD
C OD
C
COMMENT 520: CHOOSE ITERATION OPTION
C
IF( MTYPE.LT.2 ) GO TO 540
GOTO 560
C THEN
C
COMMENT 540: FIXED POINT ITERATION
C
540 DO 550 J=1 , N
TEMP=REALS(FSPACE+N+J)
REALS(YSPACE+N+J)=TEMP
IF(M.EQ.1)GO TO 550
DO 545 I=2 , M
REALS(YSPACE+N*I+J)=
+ REALS(FSPACE+N*I+J)-REALS(FSPACE+N*I-N+J)
545 TEMP=TEMP + REALS(FSPACE+N*I+J)
C OD
550 REALS(YSPACE+N*M+N+J)=REALS(FSPACE+N*M+N+J) - TEMP
C OD
GO TO 660
C ELSE
C
COMMENT 560: NEWTON ITERATION
C
560 REL= 2.E0/(1.E0+H/HFACTS)
MP1 = M+1
DO 640 I=1 , MP1
DO 570 J=1 , N
K=INTGRS(IPERM+J)
TEMP= REALS(YSPACE+N*I+K)-REALS(FSPACE+N*I+K)
IF(I .NE. 1) TEMP = TEMP+REALS(YLAST+K)
IF(J.EQ.1)GO TO 570
JM1 = J-1
DO 565 K=1 , JM1
565 TEMP=TEMP - REALS(FACTS+N*K+J) * REALS(YTEMP+K)
C OD
570 REALS(YTEMP+J)= TEMP/REALS(FACTS+N*J+J)
C OD
DO 640 NP1MJ = 1 , N
J = N + 1 -NP1MJ
TEMP=REALS(YTEMP+J)
IF(J.EQ.N)GO TO 580
JP1 = J + 1
DO 575 K = JP1 , N
575 TEMP = TEMP - REALS(FACTS+N*K+J) * REALS(YTEMP+K)
C OD
580 REALS(YTEMP+J)=TEMP
IF(I.LE.M) GO TO 585
GO TO 600
C THEN
585 IF(I.EQ.1) GO TO 590
GO TO 595
C THEN
590 REALS(FTOTAL+J)=REALS(FSPACE+N+J)
GO TO 600
C ELSE
595 REALS(FTOTAL+J)=REALS(FSPACE+N*I+J)+
+ REALS(FTOTAL+J)
C FI
C FI
600 IF(I.EQ.1 .OR. I.EQ.M+1) GOTO 601
GOTO 605
C THEN
601 REALS(W+J) = TEMP
GOTO 610
C ELSE
605 REALS(W+J)=TEMP+REALS(W+J)
C FI
610 IF(I.LE.M) GOTO 611
GOTO 625
C THEN
611 IF(I.EQ.M) GOTO 615
GOTO 620
C THEN
615 REALS(YLAST+J)=REALS(FTOTAL+J)
GOTO 625
C ELSE
620 REALS(YLAST+J)=REALS(FSPACE+N*I+J)
+ -REALS(W+J)
C FI
C FI
625 REALS(YSPACE+N*I+J)=REALS(YSPACE+N*I+J)
+ -REL*REALS(W+J)
IF(I.LE.M) GOTO 630
GOTO 640
C THEN
630 IF(I.EQ.1) GOTO631
GOTO 635
C THEN
631 REALS(FSPACE+N*I+J)=REALS(YSPACE+N*I+J)
GOTO 640
C ELSE
635 REALS(FSPACE+N*I+J)=REALS(YSPACE+N*I+J)
+ +REALS(FSPACE+N*I-N+J)
C FI
640 CONTINUE
C FI
C OD
C OD
C FI 540
C
COMMENT 660: RETRANSFORM USING T
C
660 DO 675 J=1 , N
DO 670 I=1 , M
TEMP = REALS(YSPACE+N+J)
IF(M.EQ.1)GO TO 670
DO 665 K=2 , M
665 TEMP= TEMP + REALS(T+MAXM*K+I)*REALS(YSPACE+N*K+J)
670 REALS(YVECT+I)= TEMP
C OD
C OD
DO 675 I=1 , M
675 REALS(YSPACE+N*I+J)=REALS(YVECT+I)
C OD
C OD
C
COMMENT 680: TEST CONVERGENCE
C
ITLIM = MAXITS + M - MLAST
IF( NITS.GT.1 ) XNORM=ENORM
ENORM=0.E0
DO 685 J=1 , N
TEMP=
+ ABS(REALS(YSPACE+N*ISAVE+J)-REALS(YSAVE+J))/EPS(J)
685 IF( TEMP .GT. ENORM ) ENORM=TEMP
C OD
IF( MERPST.EQ.0 )
+ ENORM=ENORM/(ABS(H)*REALS(XI+MROOTS+ISAVE))
IF(MTYPE.EQ.0 .OR. NITS.EQ.1) GOTO 690
GOTO 695
C THEN
690 DIVGD=.FALSE.
GOTO 700
C ELSE
695 DIVGD=ENORM .GT. XNORM + 1.0E0
C FI
700 IF(DIVGD) GOTO 701
GOTO 705
C THEN
701 CONVGD=.FALSE.
EXCESS=.FALSE.
GOTO 725
C ELSE
705 IF(ITLIM.EQ.NITS .AND. MTYPE.EQ.0) GOTO 706
GOTO 710
C THEN
706 CONVGD=.TRUE.
GOTO 720
C ELSE
710 IF(NITS.EQ.1) GOTO 711
GOTO 715
C THEN
711 CONVGD=ENORM.LT.1.E0
GOTO 720
C ELSE
715 CONVGD=ENORM**2.LT.XNORM.OR.ENORM.LT.1.E0
C FI
C FI
720 IF(.NOT.CONVGD) EXCESS=ITLIM.EQ.NITS
C FI
725 GO TO 465
C FI 480
C
COMMENT 740: CONVERGENCE FAILURE: TAKE REMEDIAL ACTION
C
740 IF( .NOT. CONVGD ) GO TO 741
GOTO 780
C THEN
741 STP=H*REALS(XI+MROOTS+ISAVE)
IF( MTYPE.LT.2.OR.X.EQ.XJAC )GO TO 745
GOTO 755
C THEN
745 IF (DIVGD) GOTO 746
GOTO 750
C THEN
746 H=0.9E0*H*XNORM**(1.E0-1.E0/ITLIM)/ENORM
GOTO 755
C ELSE
750 H=0.9E0*H*ENORM**(-1.E0/ITLIM)
C FI
C FI
755 IF( DIVGD ) MLAST=0
REVAL = MTYPE.GT.1
MSTATE = 0
INDEX = 5
IF(DIVGD)INDEX = 4
GO TO 060
C GO TO AND RETURN FROM PAUSE ROUTINE
760 GO TO 975
C ELSE
C
COMMENT 780: CONVERGENCE SUCCESS: SELECT APPROPRIATE ORDER
C
780 DO 785 J = 1 , N
785 REALS(FSPACE+N*M+N+J) = REALS(YSPACE+N*M+N+J)
C OD
STP=0.E0
DO 810 MP1MI = 1 , M
I = M + 1 - MP1MI
IF (STP .NE. 0.E0 )GO TO 815
C THEN
C ESCAPE FROM DO LOOP
C ELSE
ENORM = 0.E0
DO 790 J=1 , N
790 IF(ABS(REALS(FSPACE+N*I+N+J)).GT.EPS(J)*ENORM)
+ ENORM = ABS(REALS(FSPACE+N*I+N+J))/EPS(J)
C OD
TEMP = REALS(EFACTS+INTGRS(MAXI+I+1))+1
LBD=INTGRS(MAXI+I+1)-INTGRS(MAXI+I)
DO 805 MAXMK =1 , LBD
K = LBD + 1 - MAXMK
IF( REALS(EFACTS+INTGRS(MAXI+I)+K) .GT. TEMP
+ .OR.STP.NE.0.E0)GO TO 810
C THEN
C ESCAPE FROM LOOP
C ELSE
TEMP = REALS(EFACTS+INTGRS(MAXI+I)+K)
C = TEMP
IF( MERPST .EQ. 0 ) C = C/ABS(H*REALS(XI+MROOTS+K))
OK=I.EQ.M.AND.K .EQ. INTGRS(MAXI+I + 1)
+ -INTGRS(MAXI+I)
IF( OK )GO TO 795
GOTO 800
C THEN
795 EGROW = 1.5E0
IF(C*ENORM .LT. 1.3E0*ERREST) EGROW=C*ENORM/ERREST
IF(EGROW .GT. 1.4E0) MSTATE = 0
IF( EGROW .LT. 1.E0 ) EGROW = 1.E0
C FI
800 IF( ENORM*C.LT.1.E0 ) GO TO 801
GOTO 805
C THEN
801 STP = H*REALS(XI+MROOTS+K)
MNEW = I
IF( .NOT.OK ) MSTATE = 0
C FI
C FI
C OD
805 CONTINUE
C FI
810 MROOTS=MROOTS-I+1
C OD
815 IF(STP.EQ.0.E0) GO TO 816
GOTO 825
C THEN
816 STP=H*REALS(XI+(M*(M-1)/2)+ISAVE)
MSTATE=0
REVAL=MTYPE.GT.0
INDEX = 6
GO TO 060
C GO TO AND RETURN FROM PAUSE ROUTINE
820 GO TO 840
C ELSE
825 M=MNEW
C FI
C
COMMENT 840: INTERPOLATE OUTPUT VALUE IF REQUIRED
C
840 IF(STP.NE.0.E0) GO TO 841
GOTO 880
C THEN
841 IF((X+1.01*STP-XPRINT)*(XPRINT-X).GT.0.E0)GO TO 842
GOTO 880
C THEN
842 DO 845 I=1 , N
845 REALS(YTEMP+I)=Y(I)
D=(XPRINT-X)/H
P=D
PPRIME=0.E0
DO 850 K=1 , M
PPRIME=PPRIME-P
P=P+D*PPRIME/K
DO 850 I=1 , N
850 REALS(YTEMP+I)=REALS(YTEMP+I)-
+ PPRIME*REALS(FSPACE+N*K+I)/K
C OD
C OD
DO 855 I=1,N
TEMP = Y(I)
Y(I) = REALS(YTEMP+I)
855 REALS(YTEMP+I) = TEMP
C OD
TEMP = X
X = XPRINT
XPRINT=XPRINT+(XPRINT-XLAST)
XLAST=TEMP
INDEX=2
GO TO 060
C GO TO AND RETURN FROM PAUSE ROUTINE
860 TEMP = XLAST
XLAST = X
X = TEMP
DO 865 I=1,N
865 Y(I) = REALS(YTEMP+I)
C OD
GO TO 841
C FI
C FI
C
COMMENT 880: COMPLETE COMPUTATION
C
880 IF((XPRINT-X)*H.LT.0.E0)GO TO 881
GOTO 890
C THEN
881 STP = H*REALS(XI+(M*(M-1)/2)+INTGRS(MAXI+M+1)-
+ INTGRS(MAXI+M))
MSTATE = 0
INDEX = 7
GO TO 060
C GO TO AND RETURN FROM PAUSE ROUTINE
885 CONTINUE
C FI
890 IF(STP.NE.0.E0) GO TO 891
GOTO 920
C THEN
891 D=STP/H
P=D
PPRIME=0.E0
DO 895 K=1 , M
PPRIME=PPRIME-P
P=P+D*PPRIME/K
DO 895 I=1 , N
895 REALS(YRES+I)=
+ REALS(YRES+I) - PPRIME*REALS(FSPACE+N*K+I)/K
C OD
C OD
XRES=XRES+STP
TEMP=X
X=X+XRES
XRES=XRES-(X-TEMP)
DO 900 I=1 , N
TEMP=Y(I)
Y(I)=TEMP+REALS(YRES+I)
REALS(YRES+I)=REALS(YRES+I)-(Y(I)-TEMP)
YMAG(I)= .9E0* YMAG(I)
900 IF( ABS(Y(I)).GT.YMAG(I) ) YMAG(I) = ABS(Y(I))
C OD
INDEX=3
GO TO 060
C GO TO AND RETURN FROM PAUSE ROUTINE
915 CONTINUE
C FI
C
COMMENT 920: SPECIFY ORDER AND H FOR NEXT STEP
C
920 IF(H.LE.0.E0) D=-STPMAX
IF(H.GT.0.E0) D=STPMAX
IF(STPMAX .EQ. 0.E0) D=2.E2*H
IF(STP.NE.0.E0 .AND. ABS(D).GT.5.E0*ABS(STP))
+ D=5.E0*STP
MNEW=M
IF( MSTATE .EQ. 4) GO TO 925
GOTO 935
C THEN
925 MNEW = M+1
DO 930 J=1,N
930 REALS(FSPACE+N*(M+2)+J)=
+ (REALS(FSPACE+N*(M+2)+J)-REALS(FSPACE+N*M+N+J))*H/STP
C OD
C FI
935 MLAST = M
HLAST = H
MROOTS = MNEW*(MNEW+1) / 2
SCORE = 0.E0
DO 950 MNP1MI = 1, MNEW
I = MNEW + 1 - MNP1MI
MROOTS = MROOTS-I
XIM = REALS(XI+MROOTS+INTGRS(MAXI+I+1)-INTGRS(MAXI+I))
ENORM = 0.E0
HERR = D/XIM
DO 940 J = 1 , N
940 IF(ABS(REALS(FSPACE+N*I+N+J)).GT.EPS(J)*ENORM)
+ ENORM = ABS(REALS(FSPACE+N*I+N+J))/EPS(J)
C OD
ENORM = REALS(EFACTS+INTGRS(MAXI+I+1))*ENORM
IF( MERPST.EQ.0 ) ENORM=ENORM/(ABS(H)*XIM)
TEMP =1.4E0
IF( I .EQ. MLAST ) TEMP = 2.E0
IF(I.EQ.MLAST.AND.MSTATE.GE.2)TEMP=1.2E0*EGROW
TEMP = TEMP * ENORM
IF( TEMP.GT.(H/HERR)**(I+MERPST) )
+ HERR=H/TEMP**(1.E0/(I+MERPST))
TEMP=ABS(HERR)*XIM/(I+1)
IF( TEMP .GT. SCORE ) GO TO 945
GOTO 950
C THEN
945 IF( SCORE.GT.0.E0 ) MSTATE = 0
SCORE=TEMP
HNEW=HERR
M=I
ERREST = ENORM*(HNEW/H)**(M+MERPST)
C FI
950 CONTINUE
C OD
IF(MSTATE.NE.3.OR.(M.LT.MAXM.AND.EGROW.LT.1.2E0))
+ GO TO 955
GOTO 960
C THEN
955 MSTATE = MSTATE + 1
IF(MSTATE .EQ. 5)MSTATE=2
C FI
960 IF(MSTATE .EQ. 4) GO TO 961
GOTO 970
C THEN
961 DO 965 J=1,N
965 REALS(FSPACE+N*(M+2)+J)=REALS(FSPACE+N*M+N+J)
C OD
GO TO 975
C ELSE
970 H=HNEW
C FI
C FI 740
975 GO TO 220
C OD
END