C subroutine "argseg" to compute argument variation along a line segment C amended with "alnlin()" 17/10 1991 C Sven Ivansson 1990 SUBROUTINE ARGSEG (CDETSUB, UP1,UP2, RLNDET1,CDET1, DDSUP1, $ DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX, $ DRELGOALFAKT, ARGDIFSTORE, NPLINDIM, $ NPLIN,UPLIN,ARGLIN,ALNLIN,DDSUPLIN, RLNDET2,CDET2, $ IASFCALL,IASFAIL) C cdetsub subroutine cdetsub (px,py, rlndet,cdet) gives C the function as dexp(rlndet)*cdet C up1,up2 the end points C rlndet1,cdet1 function data for the first end point C ddsup1 suggestion for (half) initial "safe start" step, C but ddsupmax will be used if ddsup1 <= 0.0d0 C ddsupmin minimal step that is handled (strict limit), C always take ddsupmin >= 0.0d0 C ddsupmax maximal step C ddsupfaktmax maximal factor for change of "planned" steps C derivdrelmax maximum tolerable "relative difference" in the C "derivative estimation at up1" for the "safe C start", always take derivdrelmax < 1.0d0 C drelmax maximum allowed relative change of the analytic C function in each step, always take drelmax < 1.0d0 C drelgoalfakt for choosing "next step" (except during the "safe C start" where 0.25d0 is used), always take C drelgoalfakt < 1.0d0 C argdifstore lower limit in degrees for storage, always C take argdifstore < 180.0d0 C ( note that "argdifstore=0.0" implies that C "everything" is stored, this may be useful to C get a "good start for dzeroholo in zerofind" ) C nplindim declared upper dim of uplin,arglin,alnlin, C ddsuplin in calling program C nplin output stored number of points (excluding the C initial up1) C uplin(0:nplin) output stored points C arglin(0:nplin) output calculated argument data in degrees for the C points C alnlin(0:nplin) output calculated "ln(abs(..))" data for the points C (note that the storage may be sparse since C it is based on the argument variation C only, via "argdifstore") C ddsuplin(1:nplin) output obtained "least planned steps" between the C points C rlndet2,cdet2 output calculated function data for the second end C point C iasfcall input/output (OBS) performed number of function C evaluations C iasfail output 0 upon successful result C 1 for "ddsup<=ddsupmin" from the "safe start" C 2 for "ddsup<=ddsupmin" later on C 3 for insufficient nplindim C 4 for "zero found on the line" C note that "ddsup" is for the "planned" steps C while "ddsupnow" is for the "actual" step C (during the "safe start" they coincide) IMPLICIT NONE INTEGER NMISSDIM PARAMETER ( NMISSDIM = 20 ) C nmiss <= nmissdim INTEGER NMISS REAL*8 SUPMISS(nmissdim), RLNDETMISS(nmissdim) COMPLEX*16 UPMISS(nmissdim), CDETMISS(nmissdim) LOGICAL NEWUP INTEGER NPLINDIM REAL*8 RLNDET1, DDSUP1, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, ARGDIFSTORE COMPLEX*16 UP1,UP2, CDET1 INTEGER NPLIN, IASFCALL,IASFAIL REAL*8 ARGLIN(0:*), ALNLIN(0:*), DDSUPLIN(*), RLNDET2 COMPLEX*16 UPLIN(0:*), CDET2 INTEGER IARGRIKT, ISLASK REAL*8 DDSUP,DDSUPNOW, SUP,SUPOLD, SUPTOT, $ RLNDET,RLNDETOLD, ARG,ARGOLD, ALN,ALNOLD, DRELGOAL, $ DREL,DRELADV, ARGLINLAST, DDSUPINF, SLASK COMPLEX*16 UPRIKT, UP,UPOLD, CDET,CDETOLD, DERIV,DERIVLAST, $ CSLASK EXTERNAL CDETSUB IF ( DDSUPMIN.LT.0.0D0 .OR. DERIVDRELMAX.GE.1.0D0 .OR. $ DRELMAX.GE.1.0D0 .OR. DRELGOALFAKT.GE.1.0D0 .OR. $ ARGDIFSTORE.GE.180.0D0 ) STOP 416 ! the opposite would be risky IF (CDET1.EQ.(0.0D0,0.0D0)) THEN IASFAIL = 4 ! a zero on the line RETURN END IF UPRIKT = UP2 - UP1 SUPTOT = CDABS(UPRIKT) IF (SUPTOT.EQ.0.0D0) THEN NPLIN = 0 UPLIN(0) = UP1 ARGLIN(0) = DATAN2D (DIMAG(CDET1),DREAL(CDET1)) ALNLIN(0) = RLNDET1 + DLOG(CDABS(CDET1)) RLNDET2 = RLNDET1 CDET2 = CDET1 IASFAIL = 0 RETURN END IF UPRIKT = UPRIKT/SUPTOT C make a "safe start" IF (DDSUP1.LE.0.0D0) THEN DDSUP = DDSUPMAX ELSE DDSUP = DMIN1(2*DDSUP1,DDSUPMAX) END IF 1 SUP = DDSUP IF (SUP.LT.SUPTOT) THEN UP = UP1 + SUP*UPRIKT ELSE DDSUP = SUPTOT SUP = SUPTOT UP = UP2 END IF IASFCALL = IASFCALL + 1 CALL CDETSUB (DREAL(UP),DIMAG(UP), RLNDET,CDET) IF (CDET.EQ.(0.0D0,0.0D0)) THEN IASFAIL = 4 ! a zero on the line RETURN END IF C relevant for overflow/underflow is that 1.7*(10**38) "=" exp(88.2) IF ( RLNDET-RLNDET1 .GT. 88.0 ) THEN DDSUP = 0.5D0 * DDSUP GOTO 1 END IF DERIVLAST = ( DEXP(RLNDET-RLNDET1)*(CDET/CDET1) - 1 ) / DDSUP NMISS = 1 ! push the "miss stack" SUPMISS(1) = SUP UPMISS(1) = UP RLNDETMISS(1) = RLNDET CDETMISS(1) = CDET DDSUP = 0.5D0 * DDSUP IF (DDSUP.LE.DDSUPMIN) THEN IASFAIL = 1 ! a very small ddsup1 (?) RETURN END IF SUP = DDSUP UP = UP1 + SUP*UPRIKT 11 IASFCALL = IASFCALL + 1 CALL CDETSUB (DREAL(UP),DIMAG(UP), RLNDET,CDET) IF (CDET.EQ.(0.0D0,0.0D0)) THEN IASFAIL = 4 ! a zero on the line RETURN END IF CSLASK = DEXP(RLNDET-RLNDET1)*(CDET/CDET1) - 1 DREL = DMAX1(DABS(DREAL(CSLASK)),DABS(DIMAG(CSLASK))) !!! cdabs(cslask) DERIV = CSLASK/DDSUP SLASK = DMAX1(DABS(DREAL(DERIVLAST)),DABS(DIMAG(DERIVLAST))) + $ DERIVDRELMAX !!!cdabs(derivlast), +derivdrelmax to avoid 0-deriv probl CSLASK = (DERIV-DERIVLAST) / SLASK SLASK = DMAX1(DABS(DREAL(CSLASK)),DABS(DIMAG(CSLASK))) !!!cdabs(cslask) IF ( SLASK.GE.DERIVDRELMAX .OR. DREL.GE.DRELMAX ) THEN IF (NMISS.LT.nmissdim) THEN NMISS = NMISS + 1 ! push the "miss stack" SUPMISS(NMISS) = SUP UPMISS(NMISS) = UP RLNDETMISS(NMISS) = RLNDET CDETMISS(NMISS) = CDET END IF IF (DDSUP.LE.DDSUPMIN) THEN IASFAIL = 1 ! a "no or zero" derivative at up1 (?) RETURN END IF SLASK = 0.1D0*DDSUP IF (DREL.GE.DRELMAX) DDSUP = DDSUP * DRELMAX/DREL DDSUP = DMAX1 ( 0.25D0*DDSUP , SLASK , DDSUPMIN ) IF (DDSUP.EQ.0.0D0) THEN IASFAIL = 1 ! a "no or zero" derivative at up1 (?) RETURN END IF SUP = DDSUP UP = UP1 + SUP*UPRIKT DERIVLAST = DERIV GOTO 11 END IF SUPOLD = 0.0D0 UPOLD = UP1 RLNDETOLD = RLNDET1 CDETOLD = CDET1 ARGOLD = DATAN2D (DIMAG(CDET1),DREAL(CDET1)) ALNOLD = RLNDET1 + DLOG(CDABS(CDET1)) NPLIN = 0 UPLIN(0) = UP1 ARGLIN(0) = ARGOLD ALNLIN(0) = ALNOLD ARGLINLAST = ARGOLD DDSUPINF = DDSUP IARGRIKT = 0 DRELGOAL = DRELGOALFAKT * DRELMAX NEWUP = .TRUE. DDSUPNOW = DDSUP DO WHILE (.TRUE.) 111 CSLASK = DEXP(RLNDET-RLNDETOLD)*(CDET/CDETOLD) DREL = DMAX1 ( DABS(DREAL(CSLASK-1)), $ DABS(DIMAG(CSLASK-1)) ) !!! cdabs(cslask-1) IF (DREL.GE.DRELMAX) THEN IF (NEWUP) THEN IF (NMISS.LT.nmissdim) THEN NMISS = NMISS + 1 ! push the "miss stack" SUPMISS(NMISS) = SUP UPMISS(NMISS) = UP RLNDETMISS(NMISS) = RLNDET CDETMISS(NMISS) = CDET END IF ELSE NMISS = NMISS + 1 ! "reactivated" in the "miss stack" END IF DDSUP = DDSUPNOW * DRELGOAL/DREL IF (DDSUP.LT.DDSUPINF) DDSUPINF = DDSUP IF (DDSUP.LE.DDSUPMIN) THEN IASFAIL = 2 ! a zero on the line (?) RETURN END IF DDSUPNOW = DDSUP SUP = SUPOLD + DDSUPNOW UP = UP1 + SUP*UPRIKT IASFCALL = IASFCALL + 1 NEWUP = .TRUE. CALL CDETSUB (DREAL(UP),DIMAG(UP), RLNDET,CDET) IF (CDET.EQ.(0.0D0,0.0D0)) THEN IASFAIL = 4 ! a zero on the line RETURN END IF GOTO 111 END IF DRELADV = DREL / $ DMAX1(DABS(DREAL(CSLASK)),DABS(DIMAG(CSLASK))) !!! cdabs(cslask) ARG = DATAN2D (DIMAG(CDET),DREAL(CDET)) ARG = ARG + 360.0D0*IDNINT((ARGOLD-ARG)/360.0D0) ALN = RLNDET + DLOG(CDABS(CDET)) IF (IARGRIKT.NE.0) THEN IF (ARG.GT.ARGOLD) THEN ISLASK = 1 ELSE IF (ARG.LT.ARGOLD) THEN ISLASK = -1 ELSE ISLASK = IARGRIKT END IF IF (ISLASK.NE.IARGRIKT) THEN NPLIN = NPLIN + 1 IF (NPLIN.GT.NPLINDIM) THEN IASFAIL = 3 ! a zero on the line in this case too (?) RETURN END IF UPLIN(NPLIN) = UPOLD ARGLIN(NPLIN) = ARGOLD ALNLIN(NPLIN) = ALNOLD DDSUPLIN(NPLIN) = DDSUPINF ARGLINLAST = ARGOLD DDSUPINF = DDSUP IARGRIKT = 0 END IF END IF IF (IARGRIKT.EQ.0) THEN IF (ARG.GT.ARGOLD) THEN IARGRIKT = 1 ELSE IF (ARG.LT.ARGOLD) THEN IARGRIKT = -1 END IF END IF IF ( DABS(ARG-ARGLINLAST).GE.ARGDIFSTORE .OR. SUP.EQ.SUPTOT ) THEN NPLIN = NPLIN + 1 IF (NPLIN.GT.NPLINDIM) THEN IASFAIL = 3 ! a zero on the line in this case too (?) RETURN END IF UPLIN(NPLIN) = UP ARGLIN(NPLIN) = ARG ALNLIN(NPLIN) = ALN DDSUPLIN(NPLIN) = DDSUPINF ARGLINLAST = ARG DDSUPINF = DDSUPMAX IARGRIKT = 0 END IF IF (SUP.EQ.SUPTOT) THEN RLNDET2 = RLNDET CDET2 = CDET IASFAIL = 0 RETURN END IF SUPOLD = SUP UPOLD = UP RLNDETOLD = RLNDET CDETOLD = CDET ARGOLD = ARG ALNOLD = ALN IF (DRELADV.GT.0.0D0) THEN DDSUP = DMIN1 ( DDSUPMAX, DDSUPFAKTMAX*DDSUP, $ DDSUPNOW*DRELGOAL/DRELADV ) ELSE DDSUP = DMIN1 ( DDSUPMAX, DDSUPFAKTMAX*DDSUP ) END IF IF (DDSUP.LT.DDSUPINF) DDSUPINF = DDSUP IF (DDSUP.LE.DDSUPMIN) THEN IASFAIL = 2 ! a zero on the line (?) RETURN END IF SUP = SUPOLD + DDSUP IF (NMISS.EQ.0) THEN NEWUP = .TRUE. ELSE IF (SUPMISS(NMISS).LE.SUP) THEN NEWUP = .FALSE. ELSE NEWUP = .TRUE. END IF END IF IF (NEWUP) THEN IF (SUP.LT.SUPTOT) THEN DDSUPNOW = DDSUP UP = UP1 + SUP*UPRIKT ELSE DDSUPNOW = SUPTOT - SUPOLD SUP = SUPTOT UP = UP2 END IF IASFCALL = IASFCALL + 1 CALL CDETSUB (DREAL(UP),DIMAG(UP), RLNDET,CDET) IF (CDET.EQ.(0.0D0,0.0D0)) THEN IASFAIL = 4 ! a zero on the line RETURN END IF ELSE SUP = SUPMISS(NMISS) DDSUPNOW = SUP - SUPOLD UP = UPMISS(NMISS) RLNDET = RLNDETMISS(NMISS) CDET = CDETMISS(NMISS) NMISS = NMISS - 1 ! pop the "miss stack" END IF END DO END