C ALGORITHM 785, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 24,NO. 3, September, 1998, P. 317--333. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Fortran77/ # Fortran77/Drivers/ # Fortran77/Drivers/Dp/ # Fortran77/Drivers/Dp/RES # Fortran77/Drivers/Dp/data # Fortran77/Drivers/Dp/driver.f # Fortran77/Src/ # Fortran77/Src/Dp/ # Fortran77/Src/Dp/src.f # This archive created: Tue Mar 23 08:54:59 1999 export PATH; PATH=/bin:$PATH if test ! -d 'Fortran77' then mkdir 'Fortran77' fi cd 'Fortran77' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'RES' then echo shar: will not over-write existing file "'RES'" else cat << SHAR_EOF > 'RES' TYPE 1,2,3,4,5,6 OR 7 TO SELECT AN EXAMPLE AVAILABLE SELECTION: 1 A SQUARE-SYMMETRIC REGION 2 A MILDLY CROWDED INFINITE REGION 3 A HEAVILY CROWDED REGION 4 A CHINESE-CHARACTER-STRUCTURED REGION 5 A 4-DIRECTION-INFINITE REGION 6 AN EX. GIVEN BY REFEREE FOR CHECKING INV. MAP 7 THE UPPER HALF PLANE WITH A HORIZONTAL SLIT INPUT NPTQ(THE NUMBER OF G-J POINTS) (RECOMMENDED VALUES FOR NPTQ ARE: 2,3,4,5,6,7,OR 8) *** INPUTS ARE CHECKED WITH NO ERROR BEING FOUND *** SOLVE THE NONLINEAR SYSTEM? (1 FOR "YES", 2 FOR "NO") WANT SCREEN-DISPLAY OF THE RESIDUAL OF THE SYSTEM AS THE ITERATION GOES ON ? (1 FOR "YES", 2 FOR "NO") ***TIME FOR OBTAINING A CONVERGED RESULT MAY RANGE FROM SEVERAL SECONDS TO SEVERAL MINUTES OR LONGER DEPENDING ON THE COMPLEXITY OF THE GEOMETRY OF THE REGION AND THE LOCAL COMPUTING SYSTEM. SO, BE PATIENT ! INFO= 1 PARAMETERS DEFINING MAP: (M= 4) (N= 4) (NPTQ= 7) (TOL= 0.1E-09) U= 0.4559381277686 C=( 1.8134230823216,-1.8134230823303) K W0(K) PHI0(K) --- ----- ------- 1 ( 0.0000000000069, 1.0000000000000) 1.570796326788 2 (-1.0000000000000,-0.0000000000048) 3.141592653595 3 ( 0.0000000000118,-1.0000000000000) 4.712388980396 4 ( 1.0000000000000, 0.0000000000000) 0.000000000000 K W1(K) PHI1(K) --- ----- ------- 1 ( 0.3223969419470, 0.3223969419463) 0.7853981633963 2 (-0.3223969419468, 0.3223969419466) 2.3561944901927 3 (-0.3223969419448,-0.3223969419486) 3.9269908169932 4 ( 0.3223969419481,-0.3223969419452) -0.7853981633930 ACCURACY TEST: MAXIMUM ERROR= 0.191E-10 ACHIEVED AT KMAX= 2 MINIMUM ERROR= 0.100E-10 ACHIEVED AT KMAX= 3 WANT TO INVERSE THE MAP? (1 FOR "YES",2 FOR "NO") WANT TO TRY ANOTHER EXAMPLE? (1 FOR "YES",2 FOR "NO") TYPE 1,2,3,4,5,6 OR 7 TO SELECT AN EXAMPLE AVAILABLE SELECTION: 1 A SQUARE-SYMMETRIC REGION 2 A MILDLY CROWDED INFINITE REGION 3 A HEAVILY CROWDED REGION 4 A CHINESE-CHARACTER-STRUCTURED REGION 5 A 4-DIRECTION-INFINITE REGION 6 AN EX. GIVEN BY REFEREE FOR CHECKING INV. MAP 7 THE UPPER HALF PLANE WITH A HORIZONTAL SLIT INPUT NPTQ(THE NUMBER OF G-J POINTS) (RECOMMENDED VALUES FOR NPTQ ARE: 2,3,4,5,6,7,OR 8) *** INPUTS ARE CHECKED WITH NO ERROR BEING FOUND *** SOLVE THE NONLINEAR SYSTEM? (1 FOR "YES", 2 FOR "NO") WANT SCREEN-DISPLAY OF THE RESIDUAL OF THE SYSTEM AS THE ITERATION GOES ON ? (1 FOR "YES", 2 FOR "NO") ***TIME FOR OBTAINING A CONVERGED RESULT MAY RANGE FROM SEVERAL SECONDS TO SEVERAL MINUTES OR LONGER DEPENDING ON THE COMPLEXITY OF THE GEOMETRY OF THE REGION AND THE LOCAL COMPUTING SYSTEM. SO, BE PATIENT ! INFO= 1 PARAMETERS DEFINING MAP: (M=12) (N= 6) (NPTQ= 8) (TOL= 0.1E-09) U= 0.4280124243898 C=(-0.5933837786172, 0.3918657659855) K W0(K) PHI0(K) --- ----- ------- 1 (-0.6190110814085, 0.7853822515779) 2.238279248468 2 (-0.6230675595641, 0.7821680230097) 2.243454800623 3 (-0.6230993754454, 0.7821426777255) 2.243495477812 4 (-0.6282872591546, 0.7779814393570) 2.250146057123 5 (-0.9786505288556, 0.2055313659046) 2.934586019489 6 (-0.9335712511816,-0.3583918511451) 3.508137398713 7 (-0.7302223494245,-0.6832095728259) 3.893741637986 8 (-0.5577586554076,-0.8300031821131) 4.120706043113 9 (-0.5576885122132,-0.8300503137433) 4.120790550263 10 (-0.5543694995197,-0.8322706639083) 4.124783772183 11 (-0.3265970215046,-0.9451636818797) 4.379688067292 12 ( 1.0000000000000, 0.0000000000000) 0.000000000000 K W1(K) PHI1(K) --- ----- ------- 1 ( 0.2101301204192, 0.3728806349553) 1.0576233605451 2 ( 0.0717410804160, 0.4219571694056) 1.4023869100452 3 (-0.4122660002761, 0.1150277377349) 2.8694990601434 4 (-0.0812736520819,-0.4202252121295) 4.5213427457428 5 ( 0.2800464274801,-0.3236798323771) 5.4256390982226 6 ( 0.4239071523763, 0.0591384950455) 0.1386134964013 ACCURACY TEST: MAXIMUM ERROR= 0.811E-09 ACHIEVED AT KMAX= 2 MINIMUM ERROR= 0.223E-10 ACHIEVED AT KMAX= 1 WANT TO INVERSE THE MAP? (1 FOR "YES",2 FOR "NO") INPUT THE POINT INSIDE THE DOMAIN IN THE FORM OF ( X-COORD., Y-COORD. ) THE PREIMAGE OF ZZ= (-0.623005768209842,0.782172159414472) CHECK BY MAPPING THE PREIMAGE BACK THE POINT ENTERED= (0.499999998638322,0.499999999326471) WANT TO TRY ANOTHER PT ? (1 FOR "YES",2 FOR "NO") WANT TO TRY ANOTHER EXAMPLE? (1 FOR "YES",2 FOR "NO") TYPE 1,2,3,4,5,6 OR 7 TO SELECT AN EXAMPLE AVAILABLE SELECTION: 1 A SQUARE-SYMMETRIC REGION 2 A MILDLY CROWDED INFINITE REGION 3 A HEAVILY CROWDED REGION 4 A CHINESE-CHARACTER-STRUCTURED REGION 5 A 4-DIRECTION-INFINITE REGION 6 AN EX. GIVEN BY REFEREE FOR CHECKING INV. MAP 7 THE UPPER HALF PLANE WITH A HORIZONTAL SLIT INPUT NPTQ(THE NUMBER OF G-J POINTS) (RECOMMENDED VALUES FOR NPTQ ARE: 2,3,4,5,6,7,OR 8) *** INPUTS ARE CHECKED WITH NO ERROR BEING FOUND *** SOLVE THE NONLINEAR SYSTEM? (1 FOR "YES", 2 FOR "NO") WANT SCREEN-DISPLAY OF THE RESIDUAL OF THE SYSTEM AS THE ITERATION GOES ON ? (1 FOR "YES", 2 FOR "NO") ***TIME FOR OBTAINING A CONVERGED RESULT MAY RANGE FROM SEVERAL SECONDS TO SEVERAL MINUTES OR LONGER DEPENDING ON THE COMPLEXITY OF THE GEOMETRY OF THE REGION AND THE LOCAL COMPUTING SYSTEM. SO, BE PATIENT ! INFO= 1 PARAMETERS DEFINING MAP: (M=11) (N= 6) (NPTQ= 8) (TOL= 0.1E-09) U= 0.5022139028275 C=( 0.4495485170805, 0.9332530061513) K W0(K) PHI0(K) --- ----- ------- 1 ( 0.4783020485418, 0.8781953941810) 1.072076089882 2 (-0.6571395336034, 0.7537689522496) 2.287813895260 3 (-0.6674256911660, 0.7446764040652) 2.301542785614 4 (-0.7988453306875, 0.6015364807215) 2.496169557448 5 (-0.9813485723828,-0.1922367797332) 3.335033585256 6 (-0.9593804934419,-0.2821153466284) 3.427590959182 7 (-0.9591213109369,-0.2829952489118) 3.428508239760 8 (-0.9591213094392,-0.2829952539878) 3.428508245053 9 (-0.9588639952683,-0.2838658813208) 3.429416101094 10 (-0.7355651877502,-0.6774539501472) 3.885888394658 11 ( 1.0000000000000, 0.0000000000000) 0.000000000000 K W1(K) PHI1(K) --- ----- ------- 1 ( 0.2403845885364, 0.4409467697890) 1.0716800530526 2 (-0.2756891411992, 0.4197788722864) 2.1519010344455 3 (-0.3527528425814, 0.3574692102041) 2.3495538933856 4 (-0.3576966667014, 0.3525221962145) 2.3634801019986 5 (-0.4241878817812,-0.2688558073450) 3.7065043429625 6 ( 0.5022111128429, 0.0016740162215) 0.0033332795289 ACCURACY TEST: MAXIMUM ERROR= 0.948E-06 ACHIEVED AT KMAX= 8 MINIMUM ERROR= 0.130E-10 ACHIEVED AT KMAX= 2 WANT TO INVERSE THE MAP? (1 FOR "YES",2 FOR "NO") WANT TO TRY ANOTHER EXAMPLE? (1 FOR "YES",2 FOR "NO") TYPE 1,2,3,4,5,6 OR 7 TO SELECT AN EXAMPLE AVAILABLE SELECTION: 1 A SQUARE-SYMMETRIC REGION 2 A MILDLY CROWDED INFINITE REGION 3 A HEAVILY CROWDED REGION 4 A CHINESE-CHARACTER-STRUCTURED REGION 5 A 4-DIRECTION-INFINITE REGION 6 AN EX. GIVEN BY REFEREE FOR CHECKING INV. MAP 7 THE UPPER HALF PLANE WITH A HORIZONTAL SLIT INPUT NPTQ(THE NUMBER OF G-J POINTS) (RECOMMENDED VALUES FOR NPTQ ARE: 2,3,4,5,6,7,OR 8) *** INPUTS ARE CHECKED WITH NO ERROR BEING FOUND *** SOLVE THE NONLINEAR SYSTEM? (1 FOR "YES", 2 FOR "NO") WANT SCREEN-DISPLAY OF THE RESIDUAL OF THE SYSTEM AS THE ITERATION GOES ON ? (1 FOR "YES", 2 FOR "NO") ***TIME FOR OBTAINING A CONVERGED RESULT MAY RANGE FROM SEVERAL SECONDS TO SEVERAL MINUTES OR LONGER DEPENDING ON THE COMPLEXITY OF THE GEOMETRY OF THE REGION AND THE LOCAL COMPUTING SYSTEM. SO, BE PATIENT ! INFO= 1 PARAMETERS DEFINING MAP: (M= 4) (N=17) (NPTQ= 4) (TOL= 0.1E-09) U= 0.5080637834205 C=(-0.8252185030582, 0.7346741621724) K W0(K) PHI0(K) --- ----- ------- 1 ( 0.1160105274551, 0.9932479838991) 1.454523990458 2 (-0.9999980060311, 0.0019969811680) 3.139595671095 3 (-0.1134068924428,-0.9935486282747) 4.598737580188 4 ( 1.0000000000000, 0.0000000000000) 0.000000000000 K W1(K) PHI1(K) --- ----- ------- 1 ( 0.4647382191519, 0.2052978219151) 0.4159715655058 2 ( 0.4576198779787, 0.2207098894528) 0.4493872929813 3 ( 0.3795224218977, 0.3377743911259) 0.7272618582239 4 ( 0.2723085014565, 0.4289252709483) 1.0051364053220 5 ( 0.2578263039247, 0.4377835138834) 1.0385521273709 6 ( 0.1081199397957, 0.4964261139809) 1.3563486311691 7 (-0.5058226132834,-0.0476685631704) 3.2355548301266 8 (-0.4763529344161,-0.1766824549772) 3.4967697831037 9 (-0.4567900505999,-0.2224222509023) 3.5947252240323 10 (-0.4533269400240,-0.2293981113088) 3.6100545143865 11 (-0.4531551749891,-0.2297372312103) 3.6108027252558 12 (-0.4511594952261,-0.2336320138400) 3.6194164493863 13 (-0.3780577924477,-0.3394128954432) 3.8731801310193 14 (-0.2710041136013,-0.4297506002727) 4.1497671944776 15 (-0.2565200110852,-0.4385502159803) 4.1831260467214 16 (-0.1067624140812,-0.4967198354839) 4.5006750720130 17 ( 0.5056173096609, 0.0497990380914) 0.0981749265419 ACCURACY TEST: MAXIMUM ERROR= 0.172E-06 ACHIEVED AT KMAX= 3 MINIMUM ERROR= 0.539E-07 ACHIEVED AT KMAX= 1 WANT TO INVERSE THE MAP? (1 FOR "YES",2 FOR "NO") WANT TO TRY ANOTHER EXAMPLE? (1 FOR "YES",2 FOR "NO") TYPE 1,2,3,4,5,6 OR 7 TO SELECT AN EXAMPLE AVAILABLE SELECTION: 1 A SQUARE-SYMMETRIC REGION 2 A MILDLY CROWDED INFINITE REGION 3 A HEAVILY CROWDED REGION 4 A CHINESE-CHARACTER-STRUCTURED REGION 5 A 4-DIRECTION-INFINITE REGION 6 AN EX. GIVEN BY REFEREE FOR CHECKING INV. MAP 7 THE UPPER HALF PLANE WITH A HORIZONTAL SLIT INPUT NPTQ(THE NUMBER OF G-J POINTS) (RECOMMENDED VALUES FOR NPTQ ARE: 2,3,4,5,6,7,OR 8) *** INPUTS ARE CHECKED WITH NO ERROR BEING FOUND *** SOLVE THE NONLINEAR SYSTEM? (1 FOR "YES", 2 FOR "NO") WANT SCREEN-DISPLAY OF THE RESIDUAL OF THE SYSTEM AS THE ITERATION GOES ON ? (1 FOR "YES", 2 FOR "NO") ***TIME FOR OBTAINING A CONVERGED RESULT MAY RANGE FROM SEVERAL SECONDS TO SEVERAL MINUTES OR LONGER DEPENDING ON THE COMPLEXITY OF THE GEOMETRY OF THE REGION AND THE LOCAL COMPUTING SYSTEM. SO, BE PATIENT ! INFO= 1 PARAMETERS DEFINING MAP: (M=12) (N= 4) (NPTQ= 3) (TOL= 0.1E-09) U= 0.3732593282942 C=(-0.4884261138328,-0.8609207012823) K W0(K) PHI0(K) --- ----- ------- 1 ( 0.6542535856529, 0.7562752446433) 0.857601089415 2 ( 0.2115395646982, 0.9773694350486) 1.357646422167 3 ( 0.0984674961776, 0.9951402675987) 1.472169011342 4 (-0.5708551758343, 0.8210507707951) 2.178343368006 5 (-0.9227803940901, 0.3853262828863) 2.746031270415 6 (-0.9510906909506, 0.3089117958044) 2.827543998566 7 (-0.9984302522853,-0.0560092074704) 3.197631186249 8 ( 0.2649022418819,-0.9642752730657) 4.980491521406 9 ( 0.5034874374818,-0.8640025464650) 5.240019405018 10 ( 0.7670211252277,-0.6416218461481) 5.586574432807 11 ( 0.9554686968802,-0.2950924758139) 5.983633004144 12 ( 1.0000000000000, 0.0000000000000) 0.000000000000 K W1(K) PHI1(K) --- ----- ------- 1 (-0.2579461752072, 0.2697893564518) 2.3337567482261 2 (-0.2889380091518,-0.2362992869775) 3.8271054023731 3 ( 0.2825040966070,-0.2439548350801) 5.5708803774311 4 ( 0.2195123405734, 0.3018888181014) 0.9420935621688 ACCURACY TEST: MAXIMUM ERROR= 0.110E-03 ACHIEVED AT KMAX= 6 MINIMUM ERROR= 0.510E-05 ACHIEVED AT KMAX= 1 WANT TO INVERSE THE MAP? (1 FOR "YES",2 FOR "NO") WANT TO TRY ANOTHER EXAMPLE? (1 FOR "YES",2 FOR "NO") TYPE 1,2,3,4,5,6 OR 7 TO SELECT AN EXAMPLE AVAILABLE SELECTION: 1 A SQUARE-SYMMETRIC REGION 2 A MILDLY CROWDED INFINITE REGION 3 A HEAVILY CROWDED REGION 4 A CHINESE-CHARACTER-STRUCTURED REGION 5 A 4-DIRECTION-INFINITE REGION 6 AN EX. GIVEN BY REFEREE FOR CHECKING INV. MAP 7 THE UPPER HALF PLANE WITH A HORIZONTAL SLIT INPUT NPTQ(THE NUMBER OF G-J POINTS) (RECOMMENDED VALUES FOR NPTQ ARE: 2,3,4,5,6,7,OR 8) *** INPUTS ARE CHECKED WITH NO ERROR BEING FOUND *** SOLVE THE NONLINEAR SYSTEM? (1 FOR "YES", 2 FOR "NO") WANT SCREEN-DISPLAY OF THE RESIDUAL OF THE SYSTEM AS THE ITERATION GOES ON ? (1 FOR "YES", 2 FOR "NO") ***TIME FOR OBTAINING A CONVERGED RESULT MAY RANGE FROM SEVERAL SECONDS TO SEVERAL MINUTES OR LONGER DEPENDING ON THE COMPLEXITY OF THE GEOMETRY OF THE REGION AND THE LOCAL COMPUTING SYSTEM. SO, BE PATIENT ! INFO= 1 PARAMETERS DEFINING MAP: (M= 7) (N= 2) (NPTQ= 5) (TOL= 0.1E-09) U= 0.1772675790003 C=(-0.7515041412364, 1.2077109426813) K W0(K) PHI0(K) --- ----- ------- 1 ( 0.2067927722779, 0.9783847654852) 1.362500593615 2 (-0.5854410170147,-0.8107150026962) 4.086965044137 3 (-0.4545201624786,-0.8907364491815) 4.240555537868 4 (-0.4484706340671,-0.8937975667786) 4.247335464934 5 (-0.3868897862706,-0.9221259638897) 4.315132656973 6 ( 0.9923283906938,-0.1236299519657) 6.159238234974 7 ( 1.0000000000000, 0.0000000000000) 0.000000000000 K W1(K) PHI1(K) --- ----- ------- 1 (-0.0842113609299,-0.1559879522750) 4.2173655802602 2 ( 0.1027420798612, 0.1444571202482) 0.9525738702698 ACCURACY TEST: MAXIMUM ERROR= 0.550E-06 ACHIEVED AT KMAX= 6 MINIMUM ERROR= 0.330E-06 ACHIEVED AT KMAX= 1 WANT TO INVERSE THE MAP? (1 FOR "YES",2 FOR "NO") WANT TO TRY ANOTHER EXAMPLE? (1 FOR "YES",2 FOR "NO") TYPE 1,2,3,4,5,6 OR 7 TO SELECT AN EXAMPLE AVAILABLE SELECTION: 1 A SQUARE-SYMMETRIC REGION 2 A MILDLY CROWDED INFINITE REGION 3 A HEAVILY CROWDED REGION 4 A CHINESE-CHARACTER-STRUCTURED REGION 5 A 4-DIRECTION-INFINITE REGION 6 AN EX. GIVEN BY REFEREE FOR CHECKING INV. MAP 7 THE UPPER HALF PLANE WITH A HORIZONTAL SLIT INPUT NPTQ(THE NUMBER OF G-J POINTS) (RECOMMENDED VALUES FOR NPTQ ARE: 2,3,4,5,6,7,OR 8) *** INPUTS ARE CHECKED WITH NO ERROR BEING FOUND *** SOLVE THE NONLINEAR SYSTEM? (1 FOR "YES", 2 FOR "NO") WANT SCREEN-DISPLAY OF THE RESIDUAL OF THE SYSTEM AS THE ITERATION GOES ON ? (1 FOR "YES", 2 FOR "NO") ***TIME FOR OBTAINING A CONVERGED RESULT MAY RANGE FROM SEVERAL SECONDS TO SEVERAL MINUTES OR LONGER DEPENDING ON THE COMPLEXITY OF THE GEOMETRY OF THE REGION AND THE LOCAL COMPUTING SYSTEM. SO, BE PATIENT ! INFO= 1 PARAMETERS DEFINING MAP: (M= 3) (N= 2) (NPTQ= 7) (TOL= 0.1E-09) U= 0.0857957283188 C=(-2.7231050437514,-0.9071249391856) K W0(K) PHI0(K) --- ----- ------- 1 (-0.8002287256225, 0.5996949113413) 2.498472851106 2 (-0.3160469144053,-0.9487435627686) 4.390829054270 3 ( 1.0000000000000, 0.0000000000000) 0.000000000000 K W1(K) PHI1(K) --- ----- ------- 1 (-0.0271154752043,-0.0813981449555) 4.3908290542696 2 ( 0.0271154752043, 0.0813981449555) 1.2492364006799 ACCURACY TEST: MAXIMUM ERROR= 0.544E-07 ACHIEVED AT KMAX= 1 MINIMUM ERROR= 0.544E-07 ACHIEVED AT KMAX= 1 WANT TO INVERSE THE MAP? (1 FOR "YES",2 FOR "NO") WANT TO TRY ANOTHER EXAMPLE? (1 FOR "YES",2 FOR "NO") SHAR_EOF fi # end of overwriting check if test -f 'data' then echo shar: will not over-write existing file "'data'" else cat << SHAR_EOF > 'data' 1 # problem 1 7 1 2 2 1 2 # problem 2 8 1 2 1 (0.5,0.5) 2 1 3 # problem 3 8 1 2 2 1 4 # problem 4 4 1 2 2 1 5 # problem 5 3 1 2 2 1 6 # problem 6 5 1 2 2 1 7 # problem 7 7 1 2 2 2 SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << SHAR_EOF > 'driver.f' C C THIS IS THE FINAL VERSION SUBMITTED TO ACM TOMS. PROGRAM DRIVER C THIS PROGRAM DRIVES DSCPACK. IT READS FROM THE FILE DSCDATA. ON C OUTPUT,IT GIVES THE COMPUTED ACCESSORY PARAMETERS.THE ACCURACY C OF THE PARAMETERS IS ALSO TESTED. C C C .. Scalars in Common .. DOUBLE PRECISION DLAM INTEGER IU C .. C .. Arrays in Common .. DOUBLE PRECISION UARY(8),VARY(3) C .. C .. Local Scalars .. DOUBLE COMPLEX C,WW,ZZ,ZZ0 DOUBLE PRECISION EPS,PI,TOL,U INTEGER I,IGUESS,INV,IPOLY,ISHAPE,ISOLV,ITRY,K,LINEARC,M,N,NPTQ C .. C .. Local Arrays .. DOUBLE COMPLEX W0(30),W1(30),Z0(30),Z1(30) DOUBLE PRECISION ALFA0(30),ALFA1(30),PHI0(30),PHI1(30),QWORK(1660) C .. C .. External Functions .. DOUBLE COMPLEX WDSC,ZDSC EXTERNAL WDSC,ZDSC C .. C .. External Subroutines .. EXTERNAL ANGLES,CHECK,DSCDATA,DSCSOLV,DSCTEST,QINIT,THDATA C .. C .. Intrinsic Functions .. INTRINSIC ABS,ACOS C .. C .. Common blocks .. COMMON /PARAM4/UARY,VARY,DLAM,IU C .. PI = ACOS(-1.D0) C C SPECIFY THE GOEMETRY OF THE POLYGON REGION(READ FROM DATA FILE): 10 PRINT '(/)' WRITE (6,FMT=*) 'TYPE 1,2,3,4,5,6 OR 7 TO SELECT AN EXAMPLE' WRITE (6,FMT=*) 'AVAILABLE SELECTION:' WRITE (6,FMT=*) '1 A SQUARE-SYMMETRIC REGION' WRITE (6,FMT=*) '2 A MILDLY CROWDED INFINITE REGION' WRITE (6,FMT=*) '3 A HEAVILY CROWDED REGION' WRITE (6,FMT=*) '4 A CHINESE-CHARACTER-STRUCTURED REGION' WRITE (6,FMT=*) '5 A 4-DIRECTION-INFINITE REGION' WRITE (6,FMT=*) '6 AN EX. GIVEN BY REFEREE FOR CHECKING INV. MAP' WRITE (6,FMT=*) '7 THE UPPER HALF PLANE WITH A HORIZONTAL SLIT' C C INITIALISING ARRAYS: DO 20 K = 1,30 ALFA0(K) = 0.D0 ALFA1(K) = 0.D0 20 CONTINUE READ (5,FMT=*) IPOLY CALL DSCDATA(IPOLY,M,N,Z0,Z1,ALFA0,ALFA1) CALL ANGLES(N,Z1,ALFA1,1) IF (IPOLY.EQ.2 .OR. IPOLY.EQ.5 .OR. IPOLY.EQ.7) GO TO 30 CALL ANGLES(M,Z0,ALFA0,0) 30 CONTINUE C C GENERATE THE GAUSS-JACOBI WEIGHTS & NODES AND CHECK THE INPUT: WRITE (6,FMT=*) 'INPUT NPTQ(THE NUMBER OF G-J POINTS)' WRITE (6,FMT=*) + '(RECOMMENDED VALUES FOR NPTQ ARE: 2,3,4,5,6,7,OR 8)' READ (5,FMT=*) NPTQ CALL QINIT(M,N,ALFA0,ALFA1,NPTQ,QWORK) ISHAPE = 0 IF (IPOLY.EQ.2 .OR. IPOLY.EQ.5 .OR. IPOLY.EQ.7) ISHAPE = 1 CALL CHECK(ALFA0,ALFA1,M,N,ISHAPE) C C SPECIFY SOME PARAMETERS OF THE CALLING SEQUENCE OF DSCSOLV: IGUESS = 1 LINEARC = 1 TOL = 1.D-10 C C SOLVE THE ACCESSORY PARAMETER PROBLEM: PRINT '(/)' WRITE (6,FMT=*) + 'SOLVE THE NONLINEAR SYSTEM? (1 FOR "YES", 2 FOR "NO")' READ (5,FMT=*) ISOLV IF (ISOLV.EQ.1) CALL DSCSOLV(TOL,IGUESS,M,N,U,C,W0,W1,PHI0,PHI1, + Z0,Z1,ALFA0,ALFA1,NPTQ,QWORK,ISHAPE, + LINEARC) C C OUTPUT WILL BE ARRANGED IN DSCSOLV. C C COMPUTE DATA FOR THETA-FUNCTION AND TEST THE ACCURACY: CALL THDATA(U) CALL DSCTEST(M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,NPTQ,QWORK) C C FOLLOWING PARAGRAPH IS FOR TESTING INVERSE EVALUATIONS: WRITE (6,FMT=*) + 'WANT TO INVERSE THE MAP? (1 FOR "YES",2 FOR "NO")' READ (5,FMT=*) INV IF (INV.EQ.1) THEN DO 50 I = 1,10 WRITE (6,FMT=*) + 'INPUT THE POINT INSIDE THE DOMAIN IN THE FORM OF' WRITE (6,FMT=*) ' ( X-COORD., Y-COORD. )' READ (5,FMT=*) ZZ EPS = 1.D-6 WW = WDSC(ZZ,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0,PHI1, + NPTQ,QWORK,EPS,1) WRITE (6,FMT=*) 'THE PREIMAGE OF ZZ=',WW IF (ABS(WW).LE.1.D-12) GO TO 40 WRITE (6,FMT=*) 'CHECK BY MAPPING THE PREIMAGE BACK' ZZ0 = ZDSC(WW,0,2,M,N,U,C,W0,W1,Z0,Z1,ALFA0,ALFA1,PHI0, + PHI1,NPTQ,QWORK,1) WRITE (6,FMT=*) 'THE POINT ENTERED=',ZZ0 40 WRITE (6,FMT=*) + 'WANT TO TRY ANOTHER PT ? (1 FOR "YES",2 FOR "NO")' READ (5,FMT=*) INV IF (INV.EQ.2) GO TO 60 50 CONTINUE END IF C PRINT '(/)' 60 WRITE (6,FMT=*) + 'WANT TO TRY ANOTHER EXAMPLE? (1 FOR "YES",2 FOR "NO")' READ (5,FMT=*) ITRY IF (ITRY.EQ.1) GO TO 10 STOP END C..................................................................... C THE FOLLOWING SUBROUTINE GENERATES DATA. C SUBROUTINE DSCDATA(IPOLY,M,N,Z0,Z1,ALFA0,ALFA1) C C SPECIFY INPUT DATA: C C CALL KILLUN C C .. Scalar Arguments .. INTEGER IPOLY,M,N C .. C .. Array Arguments .. DOUBLE COMPLEX Z0(30),Z1(30) DOUBLE PRECISION ALFA0(30),ALFA1(30) C .. C .. Local Scalars .. DOUBLE PRECISION PI,Q,S C .. C .. Intrinsic Functions .. INTRINSIC ACOS,DCMPLX,SQRT C .. PI = ACOS(-1.D0) C IF (IPOLY.EQ.1) THEN M = 4 N = 4 Q = SQRT(2.D0) Z0(1) = DCMPLX(1.D0+Q,1.D0+Q) Z0(2) = DCMPLX(-1.D0-Q,1.D0+Q) Z0(3) = DCMPLX(-1.D0-Q,-1.D0-Q) Z0(4) = DCMPLX(1.D0+Q,-1.D0-Q) Z1(1) = DCMPLX(Q,0.D0) Z1(2) = DCMPLX(0.D0,Q) Z1(3) = DCMPLX(-Q,0.D0) Z1(4) = DCMPLX(0.D0,-Q) ELSE IF (IPOLY.EQ.2) THEN M = 12 N = 6 Z0(1) = (0.6875D0,0.875D0) Z0(2) = (0.D0,0.D0) Z0(3) = (1.D0,0.D0) Z0(4) = (0.875D0,0.875D0) Z0(5) = (1.125D0,1.375D0) Z0(6) = (2,1.375D0) Z0(7) = (1.25D0,2.D0) Z0(8) = (2.25D0,2.D0) Z0(9) = (2.375D0,2.75D0) Z0(10) = (1.625D0,2.25D0) Z0(11) = (1.125D0,2.625D0) Z0(12) = (-0.5D0,2.75D0) Z1(1) = (0.375D0,1.875D0) Z1(2) = (0.5D0,2.D0) Z1(3) = (1.D0,1.5D0) Z1(4) = (0.5D0,2.1875D0) Z1(5) = (0.5D0,2.5D0) Z1(6) = Z1(4) C ALFA0(1) = 1.39169261159339475D0 ALFA0(2) = 0.28801540784794967D0 ALFA0(3) = 0.454832764699133488D0 ALFA0(4) = 1.19275085295129979D0 ALFA0(5) = 1.35241638234956651D0 ALFA0(6) = 0.D0 ALFA0(7) = 2.D0 ALFA0(8) = 0.552568456711253445D0 ALFA0(9) = 0.260264501477747753D0 ALFA0(10) = 1.39199980651013222D0 ALFA0(11) = 0.819604487273064009D0 ALFA0(12) = 0.295854728586457991D0 ELSE IF (IPOLY.EQ.3) THEN M = 11 N = 6 Z0(1) = (0.5D0,2.5D0) Z0(2) = (0.5D0,0.5D0) Z0(3) = (1.D0,0.5D0) Z0(4) = (1.D0,1.D0) Z0(5) = (1.D0,0.5D0) Z0(6) = (0.5D0,0.5D0) Z0(7) = (0.5D0,2.5D0) Z0(8) = (0.D0,2.5D0) Z0(9) = (0.D0,0.D0) Z0(10) = (2.D0,0.D0) Z0(11) = (2.D0,2.5D0) Z1(1) = (1.D0,2.D0) Z1(2) = (1.D0,1.5D0) Z1(3) = (1.D0,2.D0) Z1(4) = (1.5D0,2.D0) Z1(5) = (1.5D0,0.5D0) Z1(6) = (1.5D0,2.D0) ALFA0(1) = 1.D0/2.D0 ALFA0(2) = 1.D0/2.D0 ALFA0(3) = 1.D0/2.D0 ALFA0(4) = 2.D0 ALFA0(5) = 3.D0/2.D0 ALFA0(6) = 3.D0/2.D0 ALFA0(7) = 1.D0/2.D0 ALFA0(8) = 1.D0/2.D0 ALFA0(9) = 1.D0/2.D0 ALFA0(10) = 1.D0/2.D0 ALFA0(11) = 1.D0/2.D0 ALFA1(1) = 3.D0/2.D0 ALFA1(2) = 2.D0 ALFA1(3) = 1.D0/2.D0 ALFA1(4) = 1.D0/2.D0 ALFA1(5) = 2.D0 ALFA1(6) = 3.D0/2.D0 ELSE IF (IPOLY.EQ.4) THEN M = 4 N = 17 Z0(1) = (-1.D0,-1.D0) Z0(2) = (1.D0,-1.D0) Z0(3) = (1.D0,1.D0) Z0(4) = (-1.D0,1.D0) Z1(1) = (0.D0,0.5D0) Z1(2) = (0.D0,0.D0) Z1(3) = (-0.5D0,0.D0) Z1(4) = Z1(2) Z1(5) = (0.D0,-0.5D0) Z1(6) = (-0.5D0,-0.5D0) Z1(7) = (0.5D0,-0.5D0) Z1(8) = (0.25D0,-0.5D0) Z1(9) = (0.25D0,-0.25D0) Z1(10) = Z1(8) Z1(11) = Z1(5) Z1(12) = Z1(2) Z1(13) = (0.5D0,0.D0) Z1(14) = Z1(2) Z1(15) = Z1(1) Z1(16) = (0.5D0,0.5D0) Z1(17) = (-0.5D0,0.5D0) ELSE IF (IPOLY.EQ.5) THEN M = 12 N = 4 S = 1.D0/6.D0 Z0(1) = DCMPLX(0.D0,-8.D0*S) Z0(2) = (0.D0,0.D0) Z0(3) = (1.D0,-1.5D0) Z0(4) = (1.D0,-0.5D0) Z0(5) = (0.D0,0.D0) Z0(6) = DCMPLX(0.D0,2.D0*S) Z0(7) = (0.D0,0.D0) Z0(8) = (-1.D0,0.D0) Z0(9) = (0.D0,0.D0) Z0(10) = DCMPLX(-7.D0*S,-1.D0) Z0(11) = (0.D0,0.D0) Z0(12) = DCMPLX(-2.D0*S,-10.D0*S) ALFA0(1) = 1.75D0 ALFA0(2) = 0.D0 ALFA0(3) = 1.D0 ALFA0(4) = 1.D0 ALFA0(5) = (PI-3.5D0)/PI ALFA0(6) = 3.5D0/PI ALFA0(7) = 1.5D0 ALFA0(8) = 1.D00 ALFA0(9) = 0.D0 ALFA0(10) = 1.75D0 ALFA0(11) = 0.D0 ALFA0(12) = 1.D0 Z1(1) = DCMPLX(2.D0*S,-0.5D0) Z1(2) = DCMPLX(-S,-2.D0*S) Z1(3) = DCMPLX(-4.D0*S,-5.D0*S) Z1(4) = (0.D0,-1.D0) ELSE IF (IPOLY.EQ.6) THEN M = 7 N = 2 Z0(1) = (-2.D0,-1.D0) Z0(2) = (2.D0,-1.D0) Z0(3) = (2.D0,2.D0) Z0(4) = (-0.8D0,2.D0) Z0(5) = (1.D0,0.5D0) Z0(6) = (-1.D0,2.D0) Z0(7) = (-2.D0,2.D0) Z1(1) = (0.D0,0.D0) Z1(2) = (-1.D0,0.D0) ELSE M = 3 N = 2 Z0(1) = (1.01D0,0.D0) Z0(2) = (100.D0,100.D0) Z0(3) = (-1.01D0,0.D0) Z1(1) = (0.D0,2.D0) Z1(2) = (0.D0,1.D0) ALFA0(1) = 1.D0 ALFA0(2) = -1.D0 ALFA0(3) = 1.D0 END IF RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << SHAR_EOF > 'src.f' C C ************* DSCPACK ************************************************* * THIS IS A COLLECTION OF SUBROUTINES TO SOLVE PARAMETER PROBLEM OF * * THE SCHWARZ-CHRISTOFFEL TRANSFORMATION FOR DOUBLY CONNECTED REGIONS * * A DRIVER IS NEEDED FOR A/AN PARTICULAR PROBLEM OR APPLICATION. * * AUTHOR: * * CHENGLIE HU APRIL,1994 (REVISED JULY,1995) AT WICHITA STATE UNIV. * * A FURTHER REVISION WAS DONE AT FHSU (JULY, 1997) * * REFERENCES:1. H.DAEPPEN, DIE SCHWARZ-CRISTOFFEL ABBILDUNG FUER * * ZWEIFACH ZUSAMMENHAENGENDE GEBIETE MIT ANWENDUNGEN. * * PH.D. DISSERTATION, ETH ZUERICH. * * 2. L.N.TREFETHEN, SCPACK USER'S GUIDE(MIT REPORT 1989) * * 3. HENRICI,APPLIED & COMPUTATIONAL COMPLEX ANALYSIS,VOL.3* * 4. C.HU, APPLICATION OF COMPUTATIONAL COMPLEX ANALYSIS TO* * SOME FREE BOUNDARY AND VORTEX FLOWS.PH.D. DISS. 1995 * *********************************************************************** C C ---------------------- SUBROUTINE THDATA(U) C ---------------------- C GENERATES DATA RELATED ONLY TO INNER RADIUS C U AND USED IN COMPUTING THE THETA-FUNCTION. C C C .. Scalar Arguments .. DOUBLE PRECISION U C .. C .. Scalars in Common .. DOUBLE PRECISION DLAM INTEGER IU C .. C .. Arrays in Common .. DOUBLE PRECISION UARY(8),VARY(3) C .. C .. Local Scalars .. DOUBLE PRECISION PI INTEGER K,N C .. C .. Intrinsic Functions .. INTRINSIC ACOS,EXP,LOG C .. C .. Common blocks .. COMMON /PARAM4/UARY,VARY,DLAM,IU C .. PI = ACOS(-1.D0) IF (U.GE.0.63D0) GO TO 20 IF (U.LT.0.06D0) THEN IU = 3 ELSE IF (U.LT.0.19D0) THEN IU = 4 ELSE IF (U.LT.0.33D0) THEN IU = 5 ELSE IF (U.LT.0.45D0) THEN IU = 6 ELSE IF (U.LT.0.55D0) THEN IU = 7 ELSE IU = 8 END IF DO 10 K = 1,IU N = K**2 UARY(K) = U**N 10 CONTINUE RETURN 20 VARY(1) = EXP(PI**2/LOG(U)) DLAM = -LOG(U)/PI RETURN END DOUBLE COMPLEX C ------------------------- + FUNCTION WTHETA(U,W) C ------------------------- C EVALUATES THETA-FUNCTION AT W,WHERE U IS THE C INNER RADIUS OF THE ANNULUS.THE DEFINITION OF C THETA-FUNCTION CAN BE FOUND IN REFERENCE 3. C C C .. Scalar Arguments .. DOUBLE COMPLEX W DOUBLE PRECISION U C .. C .. Scalars in Common .. DOUBLE PRECISION DLAM INTEGER IU C .. C .. Arrays in Common .. DOUBLE PRECISION UARY(8),VARY(3) C .. C .. Local Scalars .. DOUBLE COMPLEX WT,WWN,WWN0,WWP,WWP0,ZI DOUBLE PRECISION PI INTEGER K C .. C .. Intrinsic Functions .. INTRINSIC ACOS,EXP,LOG,SQRT C .. C .. Common blocks .. COMMON /PARAM4/UARY,VARY,DLAM,IU C .. WWP = (1.D0,0.D0) WWN = (1.D0,0.D0) WTHETA = (1.D0,0.D0) IF (U.GE.0.63D0) GO TO 20 WWP0 = -W WWN0 = -1.D0/W DO 10 K = 1,IU WWP = WWP*WWP0 WWN = WWN*WWN0 WTHETA = WTHETA + UARY(K)* (WWP+WWN) 10 CONTINUE RETURN 20 PI = ACOS(-1.D0) ZI = (0.D0,1.D0) WT = -ZI*LOG(-W) IF (U.GE.0.94D0) GO TO 30 WWP = EXP(WT/DLAM) WWN = 1.D0/WWP WTHETA = WTHETA + VARY(1)* (WWP+WWN) 30 WTHETA = EXP(-WT**2/ (4.D0*PI*DLAM))*WTHETA/SQRT(DLAM) RETURN END DOUBLE COMPLEX C -------------------------------------------- + FUNCTION WPROD(W,M,N,U,W0,W1,ALFA0,ALFA1) C -------------------------------------------- C COMPUTES THE PRODUCT (D-SC INTEGRAND): C M N CPROD THETA(W/U*W0(K))**(ALFA0(K)-1)*PROD THETA(U*W/W1(K))**(ALFA1(K)-1) C K=1 K=1 C THE CALLING SEQUENCE IS EXPLAINED IN THE ABOVE FORMULA. C C C .. Scalar Arguments .. DOUBLE COMPLEX W DOUBLE PRECISION U INTEGER M,N C .. C .. Array Arguments .. DOUBLE COMPLEX W0(M),W1(N) DOUBLE PRECISION ALFA0(M),ALFA1(N) C .. C .. Scalars in Common .. DOUBLE PRECISION DLAM INTEGER IU C .. C .. Arrays in Common .. DOUBLE PRECISION UARY(8),VARY(3) C .. C .. Local Scalars .. DOUBLE COMPLEX WSUM,WTH INTEGER K C .. C .. External Functions .. DOUBLE COMPLEX WTHETA EXTERNAL WTHETA C .. C .. Intrinsic Functions .. INTRINSIC EXP,LOG C .. C .. Common blocks .. COMMON /PARAM4/UARY,VARY,DLAM,IU C .. WSUM = (0.D0,0.D0) DO 10 K = 1,M WTH = LOG(WTHETA(U,W/ (U*W0(K)))) WSUM = WSUM + (ALFA0(K)-1.D0)*WTH 10 CONTINUE DO 20 K = 1,N WTH = LOG(WTHETA(U,U*W/W1(K))) WSUM = WSUM + (ALFA1(K)-1.D0)*WTH 20 CONTINUE WPROD = EXP(WSUM) RETURN END C ---------------------------------------------- SUBROUTINE QINIT(M,N,ALFA0,ALFA1,NPTQ,QWORK) C ---------------------------------------------- C COMPUTES THE GAUSS-JACOBI NODES & WEIGHTS FOR GAUSS-JACOBI C QUADRATURE. WORK ARRAY QWORK MUST BE DIMENSIONED AT LEAST C NPTQ*(2(M+N)+3). IT IS DIVIDED UP INTO 2(M+N)+3 VECTORS OF C LENGTH NPTQ: THE FIRST M+N+1 CONTAIN QUADRATURE NODES ON C OUTPUT, THE NEXT M+N+1 CONTAIN THE CORRESPONDING WEIGHTS ON C OUTPUT,AND THE LAST ONE IS A SCRATCH VECTOR USED BY GAUSSJ. C NPTQ IS THE NUMBER OF G-J NODES (SAME AS WEIGHTS) USED. SEE C COMMENT ON ROUTINE WPROD FOR THE REST OF CALLING SEQUENCE. C C C FOR EACH FINITE VERTEX,COMPUTE NODES & WEIGHTS C FOR ONE-SIDED GAUSS-JACOBI QUADRATURE: C .. Scalar Arguments .. INTEGER M,N,NPTQ C .. C .. Array Arguments .. DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(1660) C .. C .. Local Scalars .. DOUBLE PRECISION ALPHA INTEGER INODES,ISCR,IWTS,J,K C .. C .. External Subroutines .. EXTERNAL GAUSSJ C .. ISCR = NPTQ* (2* (M+N)+2) + 1 DO 30 K = 1,M + N INODES = NPTQ* (K-1) + 1 IWTS = NPTQ* (M+N+K) + 1 IF (K.LE.M) THEN ALPHA = ALFA0(K) - 1.D0 IF (ALFA0(K).GT.0.D0) THEN CALL GAUSSJ(NPTQ,0.D0,ALPHA,QWORK(ISCR),QWORK(INODES), + QWORK(IWTS)) ELSE DO 10 J = 1,NPTQ QWORK(IWTS+J-1) = 0.D0 QWORK(INODES+J-1) = 0.D0 10 CONTINUE END IF ELSE ALPHA = ALFA1(K-M) - 1.D0 CALL GAUSSJ(NPTQ,0.D0,ALPHA,QWORK(ISCR),QWORK(INODES), + QWORK(IWTS)) END IF C C TAKE SINGULARITIES INTO ACCOUNT IN ADVANCE FOR THE C PURPOSE OF SAVING CERTAIN AMOUNT OF CALCULATION IN WQSUM: DO 20 J = 1,NPTQ QWORK(IWTS+J-1) = QWORK(IWTS+J-1)* + (1.D0+QWORK(INODES+J-1))** (-ALPHA) 20 CONTINUE 30 CONTINUE C C COMPUTE NODES & WEIGHTS FOR PURE GAUSSIAN QUADRATURE: INODES = NPTQ* (M+N) + 1 IWTS = NPTQ* (2* (M+N)+1) + 1 CALL GAUSSJ(NPTQ,0.D0,0.D0,QWORK(ISCR),QWORK(INODES),QWORK(IWTS)) RETURN END DOUBLE COMPLEX C ------------------------------------------------------------ + FUNCTION WQSUM(WA,PHIA,KWA,IC,WB,PHIB,RADIUS,M,N,U,W0,W1,ALFA0, + ALFA1,NPTQ,QWORK,LINEARC) C ------------------------------------------------------------ C CALCULATES THE COMPLEX INTEGRAL FROM WA TO WB ALONG A C LINE SEGMENT (LINEARC=0)OR A CIRCULAR ARC (LINEARC=1)WITH C POSSIBLE SINGULARITY AT WA, WHERE KWA IS THE INDEX OF WA C IN W0 (IC=0) OR IN W1 (IC=1). KWA=0 AND IC=2 ( NO OTHER C VALUES PERMITTED ) IF WA IS NOT A PREVERTEX, WHICH THEN C INDICATES THAT ONLY PURE GAUSSIAN QUARATURE IS NEEDED. C PHIA & PHIB ARE ARGUMENTS OF WA & WB RESPECTIVELY. IF C INTEGRATING ALONG A CIRCULAR ARC,RADIUS SHOULD BE ASSIGNED C TO BE EITHER 1 OR U. ANY VALUE,HOWEVER,CAN BE ASSIGNED TO C RADIUS IF INTEGRATING ALONG A LINE SEGMENT.SEE DOCUMENTATIONS C OF WPROD AND QINIT FOR THE REST OF CALLING SEQUENCE. C C C .. Scalar Arguments .. DOUBLE COMPLEX WA,WB DOUBLE PRECISION PHIA,PHIB,RADIUS,U INTEGER IC,KWA,LINEARC,M,N,NPTQ C .. C .. Array Arguments .. DOUBLE COMPLEX W0(M),W1(N) DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(NPTQ* (2* (M+N)+3)) C .. C .. Scalars in Common .. DOUBLE PRECISION DLAM INTEGER IU C .. C .. Arrays in Common .. DOUBLE PRECISION UARY(8),VARY(3) C .. C .. Local Scalars .. DOUBLE COMPLEX W,WC,WH,ZI DOUBLE PRECISION PWC,PWH INTEGER I,IOFFST,IWT1,IWT2 C .. C .. External Functions .. DOUBLE COMPLEX WPROD EXTERNAL WPROD C .. C .. Intrinsic Functions .. INTRINSIC EXP C .. C .. Common blocks .. COMMON /PARAM4/UARY,VARY,DLAM,IU C .. WQSUM = (0.D0,0.D0) C C INDEX ARRANGEMENT: IWT1 = NPTQ* (IC*M+KWA-1) + 1 IF (KWA.EQ.0) IWT1 = NPTQ* (M+N) + 1 IWT2 = IWT1 + NPTQ - 1 IOFFST = NPTQ* (M+N+1) C C COMPUTE GAUSS-JACOBI SUM(W(J)*PROD(X(J))): IF (LINEARC.EQ.1) GO TO 20 C C INTEGRATE ALONG A LINE SEGMENT: WH = (WB-WA)/2.D0 WC = (WA+WB)/2.D0 DO 10 I = IWT1,IWT2 W = WC + WH*QWORK(I) WQSUM = WQSUM + QWORK(IOFFST+I)* + WPROD(W,M,N,U,W0,W1,ALFA0,ALFA1) 10 CONTINUE WQSUM = WQSUM*WH RETURN C C INTEGRATE ALONG A CIRCULAR ARC: 20 ZI = (0.D0,1.D0) PWH = (PHIB-PHIA)/2.D0 PWC = (PHIB+PHIA)/2.D0 DO 30 I = IWT1,IWT2 W = RADIUS*EXP(ZI* (PWC+PWH*QWORK(I))) WQSUM = WQSUM + QWORK(IOFFST+I)*W* + WPROD(W,M,N,U,W0,W1,ALFA0,ALFA1) 30 CONTINUE WQSUM = WQSUM*PWH*ZI RETURN END DOUBLE COMPLEX C ------------------------------------------------------------- + FUNCTION WQUAD1(WA,PHIA,KWA,IC,WB,PHIB,RADIUS,M,N,U,W0,W1,ALFA0, + ALFA1,NPTQ,QWORK,LINEARC) C ------------------------------------------------------------- C CALCULATES THE COMPLEX INTEGRAL OF WPROD FROM WA TO WB ALONG C EITHER A CIRCULAR ARC OR A LINE-EGMENT. COMPOUND ONE-SIDED C GAUSS-JACOBI QUDRATURE IS USED.SEE SUBROUTINE WQSUM FOR THE C CALLING SEQUENCE. C C C CHECK FOR ZERO-LENGTH INTEGRAL: C .. Scalar Arguments .. DOUBLE COMPLEX WA,WB DOUBLE PRECISION PHIA,PHIB,RADIUS,U INTEGER IC,KWA,LINEARC,M,N,NPTQ C .. C .. Array Arguments .. DOUBLE COMPLEX W0(M),W1(N) DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(NPTQ* (2* (N+M)+3)) C .. C .. Scalars in Common .. DOUBLE PRECISION DLAM INTEGER IU C .. C .. Arrays in Common .. DOUBLE PRECISION UARY(8),VARY(3) C .. C .. Local Scalars .. DOUBLE COMPLEX WAA,WBB,ZI DOUBLE PRECISION PHAA,PHBB,R C .. C .. External Functions .. DOUBLE COMPLEX WQSUM DOUBLE PRECISION DIST EXTERNAL WQSUM,DIST C .. C .. Intrinsic Functions .. INTRINSIC ABS,EXP,MIN C .. C .. Common blocks .. COMMON /PARAM4/UARY,VARY,DLAM,IU C .. IF (ABS(WA-WB).GT.0.D0) GO TO 10 WQUAD1 = (0.D0,0.D0) RETURN 10 ZI = (0.D0,1.D0) IF (LINEARC.EQ.1) GO TO 30 C C LINE SEGMENT INTEGRATION PATH IS CONCERNED BELOW: C STEP1:ONE-SIDED G-J QUADRATURE FOR LEFT ENDPT WA: R = MIN(1.D0,DIST(M,N,W0,W1,WA,KWA,IC)/ABS(WB-WA)) WAA = WA + R* (WB-WA) WQUAD1 = WQSUM(WA,0.D0,KWA,IC,WAA,0.D0,0.D0,M,N,U,W0,W1,ALFA0, + ALFA1,NPTQ,QWORK,LINEARC) C C STEP2:ADJOIN INTERVALS OF PURE GAUSS QUADRATURE IF NECESSARY: 20 IF (R.EQ.1.D0) RETURN R = MIN(1.D0,DIST(M,N,W0,W1,WAA,0,IC)/ABS(WAA-WB)) WBB = WAA + R* (WB-WAA) WQUAD1 = WQUAD1 + WQSUM(WAA,0.D0,0,2,WBB,0.D0,0.D0,M,N,U,W0,W1, + ALFA0,ALFA1,NPTQ,QWORK,LINEARC) WAA = WBB GO TO 20 * * CIRCULAR ARC INTEGRATION PATH IS CONCERNED BELOW: * STEP1:ONE-SIDED G-J QUADRATURE FOR LEFT ENDPT WA: 30 R = MIN(1.D0,DIST(M,N,W0,W1,WA,KWA,IC)/ABS(WB-WA)) PHAA = PHIA + R* (PHIB-PHIA) WAA = RADIUS*EXP(ZI*PHAA) WQUAD1 = WQSUM(WA,PHIA,KWA,IC,WAA,PHAA,RADIUS,M,N,U,W0,W1,ALFA0, + ALFA1,NPTQ,QWORK,LINEARC) C C STEP2:ADJOIN INTERVALS OF PURE GAUSS QUADRATURE IF NECESSARY: 40 IF (R.EQ.1.D0) RETURN R = MIN(1.D0,DIST(M,N,W0,W1,WAA,0,IC)/ABS(WAA-WB)) PHBB = PHAA + R* (PHIB-PHAA) WBB = RADIUS*EXP(ZI*PHBB) WQUAD1 = WQUAD1 + WQSUM(WAA,PHAA,0,2,WBB,PHBB,RADIUS,M,N,U,W0,W1, + ALFA0,ALFA1,NPTQ,QWORK,LINEARC) PHAA = PHBB WAA = WBB GO TO 40 END DOUBLE COMPLEX C --------------------------------------------------------------- + FUNCTION WQUAD(WA,PHIA,KWA,ICA,WB,PHIB,KWB,ICB,RADIUS,M,N,U,W0, + W1,ALFA0,ALFA1,NPTQ,QWORK,LINEARC,IEVL) C --------------------------------------------------------------- C CALCULATES THE COMPLEX INTEGRAL OF WPROD FROM WA TO WB C ALONG A CIRCULAR ARC OR A LINE-SEGMENT.FUNCTION WQUAD1 C IS CALLED FOUR TIMES,ONE FOR EACH 1/4 OF THE INTERVAL. C NOTE: WQUAD1 ALLOWS ONLY THE LEFT ENDPOINT TO BE A C POSSIBLE SINGULARITY. SEE ROUTINE WQSUM FOR THE CALLING C SEQUENCE. C C C .. Scalar Arguments .. DOUBLE COMPLEX WA,WB DOUBLE PRECISION PHIA,PHIB,RADIUS,U INTEGER ICA,ICB,IEVL,KWA,KWB,LINEARC,M,N,NPTQ C .. C .. Array Arguments .. DOUBLE COMPLEX W0(M),W1(N) DOUBLE PRECISION ALFA0(M),ALFA1(N),QWORK(NPTQ* (2* (M+N)+3)) C .. C .. Scalars in Common .. DOUBLE PRECISION DLAM INTEGER IU C .. C .. Arrays in Common .. DOUBLE PRECISION UARY(8),VARY(3) C .. C .. Local Scalars .. DOUBLE COMPLEX WMID,WMIDA,WMIDB,WQA,WQB,ZI DOUBLE PRECISION PHMID,PHMIDA,PHMIDB,PI C .. C .. External Functions .. DOUBLE COMPLEX WQUAD1 EXTERNAL WQUAD1 C .. C .. Intrinsic Functions .. INTRINSIC ACOS,EXP C .. C .. Common blocks .. COMMON /PARAM4/UARY,VARY,DLAM,IU C .. PI = ACOS(-1.D0) ZI = (0.D0,1.D0) C C DETERMINE MIDPTS ON A LINE SEGMENT OR ON A CIRCULAR ARC: IF (LINEARC.EQ.0) THEN WMID = (WA+WB)/2.D0 WMIDA = (WA+WMID)/2.D0 WMIDB = (WB+WMID)/2.D0 PHMID = 0.D0 PHMIDA = 0.D0 PHMIDB = 0.D0 ELSE C IF (IEVL.EQ.1) GO TO 10 C IF (PHIB.LT.PHIA) PHIA = PHIA - 2.D0*PI 10 PHMID = (PHIA+PHIB)/2.D0 WMID = RADIUS*EXP(ZI*PHMID) PHMIDA = (PHIA+PHMID)/2.D0 WMIDA = RADIUS*EXP(ZI*PHMIDA) PHMIDB = (PHIB+PHMID)/2.D0 WMIDB = RADIUS*EXP(ZI*PHMIDB) END IF C C COMPOUND GAUSS-JACOBI PROCESS ACCORDING TO ONE-QUATER RULE: WQA = WQUAD1(WA,PHIA,KWA,ICA,WMIDA,PHMIDA,RADIUS,M,N,U,W0,W1, + ALFA0,ALFA1,NPTQ,QWORK,LINEARC) - + WQUAD1(WMID,PHMID,0,2,WMIDA,PHMIDA,RADIUS,M,N,U,W0,W1,ALFA0, + ALFA1,NPTQ,QWORK,LINEARC) WQB = WQUAD1(WB,PHIB,KWB,ICB,WMIDB,PHMIDB,RADIUS,M,N,U,W0,W1, + ALFA0,ALFA1,NPTQ,QWORK,LINEARC) - + WQUAD1(WMID,PHMID,0,2,WMIDB,PHMIDB,RADIUS,M,N,U,W0,W1,ALFA0, + ALFA1,NPTQ,QWORK,LINEARC) WQUAD = WQA - WQB RETURN END C ----------------------------------------------- SUBROUTINE XWTRAN(M,N,X,U,C,W0,W1,PHI0,PHI1) C ----------------------------------------------- C TRANSFORMS X(K)(UNCONSTRAINED PARAMETERS) TO ACTUAL C D-SC PARAMETERS:U,C,W0,W1.PHI0 & PHI1 ARE ARGUMENTS C OF THE PREVERTICES CONTAINED IN W0 & W1. C C C .. Scalar Arguments .. DOUBLE COMPLEX C DOUBLE PRECISION U INTEGER M,N C .. C .. Array Arguments .. DOUBLE COMPLEX W0(M),W1(N) DOUBLE PRECISION PHI0(M),PHI1(N),X(M+N+2) C .. C .. Local Scalars .. DOUBLE PRECISION DPH,PH,PHSUM,PI INTEGER I C .. C .. Intrinsic Functions .. INTRINSIC ABS,ACOS,COS,DCMPLX,EXP,SIN,SQRT C .. PI = ACOS(-1.D0) IF (ABS(X(1)).LE.1.D-14) THEN U = 0.5D0 ELSE U = (X(1)-2.D0-SQRT(0.9216D0*X(1)**2+4.D0))/ (2.D0*X(1)) U = (0.0196D0*X(1)-1.D0)/ (U*X(1)) END IF C = DCMPLX(X(2),X(3)) IF (ABS(X(N+3)).LE.1.D-14) THEN PHI1(N) = 0.D0 ELSE PH = (1.D0+SQRT(1.D0+PI*PI*X(N+3)**2))/X(N+3) PHI1(N) = PI*PI/PH END IF DPH = 1.D0 PHSUM = DPH DO 10 I = 1,N - 1 DPH = DPH/EXP(X(3+I)) PHSUM = PHSUM + DPH 10 CONTINUE DPH = 2.D0*PI/PHSUM PHI1(1) = PHI1(N) + DPH W1(1) = U*DCMPLX(COS(PHI1(1)),SIN(PHI1(1))) W1(N) = U*DCMPLX(COS(PHI1(N)),SIN(PHI1(N))) PHSUM = PHI1(1) DO 20 I = 1,N - 2 DPH = DPH/EXP(X(3+I)) PHSUM = PHSUM + DPH PHI1(I+1) = PHSUM W1(I+1) = U*DCMPLX(COS(PHSUM),SIN(PHSUM)) 20 CONTINUE DPH = 1.D0 PHSUM = DPH DO 30 I = 1,M - 1 DPH = DPH/EXP(X(N+3+I)) PHSUM = PHSUM + DPH 30 CONTINUE DPH = 2.D0*PI/PHSUM PHSUM = DPH PHI0(1) = DPH W0(1) = DCMPLX(COS(DPH),SIN(DPH)) DO 40 I = 1,M - 2 DPH = DPH/EXP(X(N+3+I)) PHSUM = PHSUM + DPH PHI0(I+1) = PHSUM W0(I+1) = DCMPLX(COS(PHSUM),SIN(PHSUM)) 40 CONTINUE RETURN END C -------------------------------------- SUBROUTINE DSCFUN(NDIM,X,FVAL,IFLAG) C -------------------------------------- C FORMS THE NONLINEAR SYSTEM SATISFIED BY D-SC PARAMETERS.THE C SUBROUTINE WILL BE CALLED BY NONLINEAR SYSTEM SOLVER HYBRD. C SEE ROUTINE DSCSOLV FOR THE VARIABLES IN THE COMMON BLOCKS. C NDIM: THE DIMENSION OF THE SYSTEM. C X: UNKNOWNS C FVAL: THE VECTOR DEFINED BY THE SYSTEM. C IFLAG:(ON OUTPUT)=1,THE ITERATION WAS SUCCESSFULLY COMPLETED. C =2,3,OR 4, UNSUCCESSFUL TERMINATION OF THE ITERATION. C =5 MAY INDICATE THE TOLERANCE GIVEN IN THE CALLING C SEQUENCE OF HYBRD IS TOO SMALL. C C C TRANSFORM UNCONSTRAINED X(K) TO ACTUAL D-SC PARAMETERS C AND COMPUTE DATA TO BE USED IN WTHETA: C C .. Scalar Arguments .. INTEGER IFLAG,NDIM C .. C .. Array Arguments .. DOUBLE PRECISION FVAL(NDIM),X(NDIM) C .. C .. Scalars in Common .. DOUBLE COMPLEX C DOUBLE PRECISION DLAM,U INTEGER ICOUNT,ISHAPE,ISPRT,IU,LINEARC,M,N,NPTQ,NSHAPE C .. C .. Arrays in Common .. DOUBLE COMPLEX W0(30),W1(30),Z0(30),Z1(30) DOUBLE PRECISION ALFA0(30),ALFA1(30),PHI0(30),PHI1(30), + QWORK(1660),UARY(8),VARY(3) INTEGER IND(20) C .. C .. Local Scalars .. DOUBLE COMPLEX WA,WAI,WARC,WB,WBI,WCIRCLE,WIN1,WIN2,WIN3,WIN4, + WINT1,WINT2,WINT3,WLINE,WLINE1,WLINE2,WX,ZI DOUBLE PRECISION FMAXN,PHIA,PHIB,RADIUS,TEST1 INTEGER I,J,K C .. C .. External Functions .. DOUBLE COMPLEX WQUAD DOUBLE PRECISION FMAX EXTERNAL WQUAD,FMAX C .. C .. External Subroutines .. EXTERNAL THDATA,XWTRAN C .. C .. Intrinsic Functions .. INTRINSIC ABS,COS,DCMPLX,DIMAG,EXP C .. C .. Common blocks .. COMMON /PARAM1/W0,W1,Z0,Z1,C COMMON /PARAM2/U,PHI0,PHI1,ALFA0,ALFA1,QWORK COMMON /PARAM3/M,N,NPTQ,ISHAPE,LINEARC,NSHAPE,IND COMMON /PARAM4/UARY,VARY,DLAM,IU COMMON /PARAM5/ISPRT,ICOUNT C .. CALL XWTRAN(M,N,X,U,C,W0,W1,PHI0,PHI1) CALL THDATA(U) ZI = (0.D0,1.D0) ICOUNT = ICOUNT + 1 C C TWO EQUATIONS TO ELIMINATE POSSIBLE ROTATION OF THE INNER POLYGON: WIN1 = Z1(1) - Z1(N) - C*WQUAD(W1(N),PHI1(N),N,1,W1(1),PHI1(1),1, + 1,U,M,N,U,W0,W1,ALFA0,ALFA1,NPTQ,QWORK,LINEARC,2) FVAL(1) = DIMAG(ZI*WIN1) FVAL(2) = DIMAG(WIN1) C C N-1 SIDE LENGTH CONDITIONS FOR THE INNER POLYGON: DO 10 I = 1,N - 1 WINT1 = WQUAD(W1(I),PHI1(I),I,1,W1(I+1),PHI1(I+1),I+1,1,U,M,N, + U,W0,W1,ALFA0,ALFA1,NPTQ,QWORK,LINEARC,2) FVAL(I+2) = ABS(Z1(I+1)-Z1(I)) - ABS(C*WINT1) 10 CONTINUE C C TWO EQUATIONS TO FIX THE RELATIVE POSITION OF THE INNER POLYGON: TEST1 = COS(PHI1(N)) IF (TEST1.GE.U) GO TO 20 C C IF THE LINE PATH FROM W0(M) TO W1(N) IS OUT OF DOMAIN,THE C COMBINATION OF TWO DIFFERENT PATHS WILL BE USED INSTEAD: WX = DCMPLX(U,0.D0) WLINE = WQUAD(W0(M),0.D0,M,0,WX,0.D0,0,2,0.D0,M,N,U,W0,W1,ALFA0, + ALFA1,NPTQ,QWORK,0,2) IF (PHI1(N).LE.0.D0) THEN WARC = WQUAD(W1(N),PHI1(N),N,1,WX,0.D0,0,2,U,M,N,U,W0,W1, + ALFA0,ALFA1,NPTQ,QWORK,1,2) WIN2 = WLINE - WARC ELSE WARC = WQUAD(WX,0.D0,0,2,W1(N),PHI1(N),N,1,U,M,N,U,W0,W1, + ALFA0,ALFA1,NPTQ,QWORK,1,2) WIN2 = WLINE + WARC END IF GO TO 30 20 WIN2 = WQUAD(W0(M),0.D0,M,0,W1(N),0.D0,N,1,0.D0,M,N,U,W0,W1,ALFA0, + ALFA1,NPTQ,QWORK,0,2) 30 FVAL(N+2) = DIMAG(ZI* (Z1(N)-Z0(M)-C*WIN2)) FVAL(N+3) = DIMAG(Z1(N)-Z0(M)-C*WIN2) C C TWO EQUATIONS TO ELIMINATE POSSIBLE ROTATION OF THE OUTER POLYGON: WIN3 = Z0(1) - Z0(M) - C*WQUAD(W0(M),PHI0(M),M,0,W0(1),PHI0(1),1, + 0,1.D0,M,N,U,W0,W1,ALFA0,ALFA1,NPTQ,QWORK,LINEARC,2) FVAL(N+4) = DIMAG(ZI*WIN3) FVAL(N+5) = DIMAG(WIN3) C IF (M.EQ.3) THEN C C CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL: IF (ISPRT.NE.1) GO TO 40 FMAXN = FMAX(NDIM,FVAL) WRITE (6,FMT=9000) ICOUNT,FMAXN 40 CONTINUE RETURN END IF IF (ISHAPE.EQ.1) GO TO 70 C C M-3 SIDE LENGTH CONDITIONS OF THE OUTER POLYGON: DO 50 J = 1,M - 3 WINT2 = WQUAD(W0(J),PHI0(J),J,0,W0(J+1),PHI0(J+1),J+1,0,1.D0, + M,N,U,W0,W1,ALFA0,ALFA1,NPTQ,QWORK,LINEARC,2) FVAL(N+5+J) = ABS(Z0(J+1)-Z0(J)) - ABS(C*WINT2) 50 CONTINUE C C CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL: IF (ISPRT.NE.1) GO TO 60 FMAXN = FMAX(NDIM,FVAL) WRITE (6,FMT=9010) ICOUNT,FMAXN 60 CONTINUE RETURN C C OUTER POLYGON CONTAINS SOME INFINITE VERTICES & FOR EACH OF THEM C TWO LENGTH CONDITIONS WILL BE REPLACED BY A COMPLEX INTEGRAL: 70 DO 100 K = 1,NSHAPE - 1 IF (IND(K+1).EQ.2 .OR. IND(K).GE.IND(K+1)-2) GO TO 90 DO 80 J = IND(K) + 1,IND(K+1) - 2 WINT3 = WQUAD(W0(J),PHI0(J),J,0,W0(J+1),PHI0(J+1),J+1,0, + 1.D0,M,N,U,W0,W1,ALFA0,ALFA1,NPTQ,QWORK,LINEARC,2) FVAL(N+5+J) = ABS(Z0(J+1)-Z0(J)) - ABS(C*WINT3) 80 CONTINUE 90 IF (K.EQ.NSHAPE-1 .OR. IND(K+1).EQ.M-1) GO TO 100 C C THE COMBINATION OF THREE DIFFERENT PATHS IS USED TO INTEGRATE C FROM WA TO WB TO AVOID DOMAIN PROBLEM.THE BY-PRODUCT OF THIS IS C THAT IT IS NUMERICALLY MORE EFFICIENT THAN A SINGLE LINE PATH: WA = W0(IND(K+1)-1) WB = W0(IND(K+1)+1) PHIA = PHI0(IND(K+1)-1) PHIB = PHI0(IND(K+1)+1) RADIUS = (1.D0+U)/2.D0 WAI = RADIUS*EXP(ZI*PHIA) WBI = RADIUS*EXP(ZI*PHIB) WLINE1 = WQUAD(WA,0.D0,IND(K+1)-1,0,WAI,0.D0,0,2,0.D0,M,N,U, + W0,W1,ALFA0,ALFA1,NPTQ,QWORK,0,2) WLINE2 = WQUAD(WB,0.D0,IND(K+1)+1,0,WBI,0.D0,0,2,0.D0,M,N,U, + W0,W1,ALFA0,ALFA1,NPTQ,QWORK,0,2) WCIRCLE = WQUAD(WAI,PHIA,0,2,WBI,PHIB,0,2,RADIUS,M,N,U,W0,W1, + ALFA0,ALFA1,NPTQ,QWORK,1,2) WIN4 = C* (WLINE1+WCIRCLE-WLINE2) FVAL(N+5+IND(K+1)-1) = DIMAG(ZI* + (Z0(IND(K+1)+1)-Z0(IND(K+1)-1)-WIN4)) FVAL(N+5+IND(K+1)) = DIMAG(Z0(IND(K+1)+1)-Z0(IND(K+1)-1)-WIN4) 100 CONTINUE C C CALCULATE THE MAXIMUM-NORM OF THE FUNCTION FVAL: IF (ISPRT.NE.1) GO TO 110 FMAXN = FMAX(NDIM,FVAL) WRITE (6,FMT=9020) ICOUNT,FMAXN 110 CONTINUE RETURN 9000 FORMAT (2X,I5,6X,E10.4) 9010 FORMAT (2X,I5,6X,E10.4) 9020 FORMAT (2X,I5,6X,E10.4) END C ----------------------------------------------------------------- SUBROUTINE DSCSOLV(TOL,IGUESS,M,N,U,C,W0,W1,PHI0,PHI1,Z0,Z1,ALFA0, + ALFA1,NPTQ,QWORK,ISHAPE,LINEARC) C ----------------------------------------------------------------- C SOLVES THE NONLINEAR SYSTEM FOR D-SC PARAMETERS. C CALLING SEQUENCE: C TOL A TOLERANCE TO CONTROL THE CONVERGENCE IN HYBRD C IGUESS (=0 )A NON-EQUALLY SPACED INITIAL GUESS OR(=1)THE C OTHER EQUALLY-SPACED INITIAL GUESS PROVIDED, OR C (=2)USER-SUPPLIED INITIAL GUESS WHICH IS REQUIRED C TO BE THE ARGUMENTS OF THE INITIAL PREVERTICES. C ROUTINE ARGUM MAY BE USED FOR COMPUTING ARGUMENTS C IF NEEDED. NOTE: C WILL BE COMPUTED IN THE ROUTINE C (NOT SUPPLIED!) C M,N,U,C,W0,W1,PHI0,PHI1,ALFA0,ALFA1,Z0,Z1 C CONSTITUTE THE GEOMETRY OF THE POLYGONAL REGION C AND THE MAPPING FUNCTION. ON RETURN U,C,W0,& W1 C WILL CONTAIN COMPUTED PARAMETERS. (PHI0 & PHI1 C WILL CONTAIN THE ARGUMENTS OF THE PREVERTICES.) C QWORK SEE CALLING SEQUENCE DOCUMENTATION IN QINIT.THE ARRAY C MUST HAVE BEEN FILLED BY QINIT BEFORE CALLING DSCSOLV. C NPTQ THE NUMBER OF GAUSS-JACOBI POINTS USED C ISHAPE INDICATES THAT OUTER POLYGON CONTAINS NO INFINITE C VERTICES (ISHAPE=0) OR IT HAS SOME INFINITE VERTICES C (ISHAPE=1). C LINEARC INTEGER VARIABLE TO CONTROL INTEGRATION PATH.IN PATICULAR: C LINEARC=0:INTEGRATING ALONG LINE SEGMENT WHENEVER POSSIBLE C LINEARC=1:INTEGRATING ALONG CIRCULAR ARC WHENEVER POSSIBLE C THE DISCRIPTION OF ARRAY IND & VARIABLE NSHAPE NEEDED IN DSCFUN: C IND CONTAINS INDICES CORRESPONDING TO THOSE INFINITE C VERTICES,BUT THE FIRST & THE LAST ENTRIS MUST BE 0 C & M-1(M IS THE # OF VERTICES ON THE OUTER POLYGON). C NSHAPE THE DIMENSION OF THE INTERGE ARRAY IND(NSHAPE) C C C .. Scalar Arguments .. DOUBLE COMPLEX C DOUBLE PRECISION TOL,U INTEGER IGUESS,ISHAPE,LINEARC,M,N,NPTQ C .. C .. Array Arguments .. DOUBLE COMPLEX W0(M),W1(N),Z0(M),Z1(N) DOUBLE PRECISION ALFA0(M),ALFA1(N),PHI0(M),PHI1(N), + QWORK(NPTQ* (2* (N+M)+3)) C .. C .. Scalars in Common .. DOUBLE COMPLEX C2 DOUBLE PRECISION DLAM,U2 INTEGER ICOUNT,ISHAPE2,ISPRT,IU,LINEARC2,M2,N2,NPTQ2,NSHAPE C .. C .. Arrays in Common .. DOUBLE COMPLEX W02(30),W12(30),Z02(30),Z12(30) DOUBLE PRECISION ALFA02(30),ALFA12(30),PHI02(30),PHI12(30), + QWORK2(1660),UARY(8),VARY(3) INTEGER IND(20) C .. C .. Local Scalars .. DOUBLE COMPLEX C1,WINT,ZI DOUBLE PRECISION AVE,BOTM,DSTEP,FACTOR,PI,TOP INTEGER I,INFO,K,KM,KN,MAXFUN,NFEV,NM,NWDIM C .. C .. Local Arrays .. DOUBLE PRECISION DIAG(42),FJAC(42,42),FVAL(42),QW(1114),X(42) C .. C .. External Functions .. DOUBLE COMPLEX WQUAD EXTERNAL WQUAD C .. C .. External Subroutines .. EXTERNAL DSCFUN,DSCPRINT,HYBRD,THDATA,XWTRAN C .. C .. Intrinsic Functions .. INTRINSIC ACOS,DBLE,DIMAG,LOG C .. C .. Common blocks .. COMMON /PARAM1/W02,W12,Z02,Z12,C2 COMMON /PARAM2/U2,PHI02,PHI12,ALFA02,ALFA12,QWORK2 COMMON /PARAM3/M2,N2,NPTQ2,ISHAPE2,LINEARC2,NSHAPE,IND COMMON /PARAM4/UARY,VARY,DLAM,IU COMMON /PARAM5/ISPRT,ICOUNT C .. ZI = (0.D0,1.D0) PI = ACOS(-1.D0) ICOUNT = 0 IF (ISHAPE.EQ.0.D0) GO TO 20 C C DETERMIN THE NUMBER OF OUTER BOUNDARY COMPONENTS,ETC.: NSHAPE = 1 DO 10 I = 2,M - 1 IF (ALFA0(I).GT.0.D0) GO TO 10 NSHAPE = NSHAPE + 1 IND(NSHAPE) = I 10 CONTINUE IND(1) = 0 NSHAPE = NSHAPE + 1 IND(NSHAPE) = M - 1 C C FIX ONE PREVERTEX: 20 W0(M) = (1.D0,0.D0) PHI0(M) = 0.D0 C C FOLLOWING TWO VALUE ASSIGNMENTS ARE TO SATISFY THE COMPILER WATFOR77: X(2) = 0.D0 X(3) = 0.D0 C INITIAL GUESS (IGUESS=0): IF (IGUESS.EQ.0) THEN X(1) = 1.D0/0.5D0 - 1.D0/0.46D0 AVE = 2.D0*PI/DBLE(N) DO 30 I = 1,N - 2 X(3+I) = LOG((AVE+0.0001D0*DBLE(I))/ + (AVE+0.0001D0*DBLE(I+1))) 30 CONTINUE X(N+2) = LOG((AVE+0.0001D0*DBLE(N-1))/ + (2.D0*PI-DBLE(N-1)* (AVE+DBLE(N)*0.00005D0))) X(N+3) = 1.D0/ (4.D0-0.1D0) - 1.D0/ (4.D0+0.1D0) AVE = 2.D0*PI/DBLE(M) DO 40 I = 1,M - 2 X(N+3+I) = LOG((AVE+0.0001D0*DBLE(I-1))/ + (AVE+0.0001D0*DBLE(I))) 40 CONTINUE X(M+N+2) = LOG((AVE+0.0001D0*DBLE(M-2))/ + (2.D0*PI-DBLE(M-1)* (AVE+DBLE(M-2)*0.00005D0))) ELSE IF (IGUESS.EQ.1) THEN C C INITIAL GUESS (IGUESS=1): X(1) = 1.D0/0.53D0 - 1.D0/0.43D0 DO 50 I = 1,N - 1 X(3+I) = 0.D0 50 CONTINUE X(N+3) = 1.D0/ (4.D0-0.1D0) - 1.D0/ (4.D0+0.1D0) DO 60 I = 1,M - 1 X(N+3+I) = 0.D0 60 CONTINUE ELSE X(1) = 1.D0/ (0.98D0-U) - 1.D0/ (U-0.02D0) DO 70 K = 1,N - 1 KN = K - 1 IF (KN.EQ.0) KN = N TOP = PHI1(K) - PHI1(KN) IF (TOP.LT.0.D0) TOP = TOP + 2.D0*PI BOTM = PHI1(K+1) - PHI1(K) IF (BOTM.LT.0.D0) BOTM = BOTM + 2.D0*PI X(3+K) = LOG(TOP) - LOG(BOTM) 70 CONTINUE X(N+3) = 1.D0/ (PI-PHI1(N)) - 1.D0/ (PI+PHI1(N)) DO 80 K = 1,M - 1 KM = K - 1 IF (KM.EQ.0) KM = M TOP = PHI0(K) - PHI0(KM) IF (TOP.LT.0.D0) TOP = TOP + 2.D0*PI BOTM = PHI0(K+1) - PHI0(K) IF (BOTM.LT.0.D0) BOTM = BOTM + 2.D0*PI X(N+3+K) = LOG(TOP) - LOG(BOTM) 80 CONTINUE END IF C C CALCULATE THE INITIAL GUESS X(2) & X(3) TO MATCH C THE CHOICE FOR X(1),X(4),...,X(M+N+2): CALL XWTRAN(M,N,X,U,C,W0,W1,PHI0,PHI1) CALL THDATA(U) WINT = WQUAD(W0(M),0.D0,M,0,W1(N),0.D0,N,1,0.D0,M,N,U,W0,W1,ALFA0, + ALFA1,NPTQ,QWORK,0,2) C1 = (Z1(N)-Z0(M))/WINT X(2) = DIMAG(ZI*C1) X(3) = DIMAG(C1) C C HYBRD CONTROL PARAMETERS: DSTEP = 1.D-7 MAXFUN = 200* (M+N) FACTOR = 2.D0 NM = M + N + 2 C C COPY RELEVANT DATA TO THE COMMON BLOCK IN DSCFUN: M2 = M N2 = N ISHAPE2 = ISHAPE LINEARC2 = LINEARC NPTQ2 = NPTQ DO 90 K = 1,M W02(K) = W0(K) PHI02(K) = PHI0(K) Z02(K) = Z0(K) ALFA02(K) = ALFA0(K) 90 CONTINUE DO 100 K = 1,N W12(K) = W1(K) PHI12(K) = PHI1(K) Z12(K) = Z1(K) ALFA12(K) = ALFA1(K) 100 CONTINUE NWDIM = NPTQ* (2* (M+N)+3) DO 110 I = 1,NWDIM QWORK2(I) = QWORK(I) 110 CONTINUE C C CHOOSE SCREEN DISPLAY: PRINT '(/)' WRITE (6,FMT=*) + 'WANT SCREEN-DISPLAY OF THE RESIDUAL OF THE SYSTEM' WRITE (6,FMT=*) + ' AS THE ITERATION GOES ON ? ' WRITE (6,FMT=*) '(1 FOR "YES", 2 FOR "NO")' PRINT '(/)' READ (5,FMT=*) ISPRT IF (ISPRT.EQ.1) THEN WRITE (6,FMT=*) + '# OF ITERATIONS>