C to find modal slownesses C note: for a rigid or free bottom the analytic function is entire C for a true half-space there are two (or only one) branch-cuts C which should be "avoided" C Sven Ivansson 1990 PROGRAM main ! RMODFAUT IMPLICIT NONE INCLUDE 'RPPARAMETER.INC' C nol <= noldim , nhyd <= nhyddim INTEGER NFRDIM, NRECTDIM, NZERODIM, ICWORKDIM, NTAB,JTAB PARAMETER ( NFRDIM = 4400 , NRECTDIM = 15 , NZERODIM = 10000 , $ ICWORKDIM = 22 , NTAB = 17 , JTAB = 9 ) C nfr <= nfrdim , nzero <= nzerodim INCLUDE 'INMODEL.INC' INCLUDE 'INPARAM.INC' INCLUDE 'DPCFSUB.INC' INCLUDE 'HARKLIQ.INC' INCLUDE 'INSPEC.INC' INCLUDE 'IADR00Z.INC' INCLUDE 'SOURCFW.INC' ! for convenience only COMMON /CORNTEST/ ICORNTEST, XYHYPA,XYHYPB INTEGER ICORNTEST REAL*8 XYHYPA,XYHYPB INTEGER MHQ0,MHQNOW,MHQEFF, ICUTTYP0,ISHEETTYP0, $ IXTYP0,NRUS0,IRUS0,JRUS0, NFR, IUROT, IRECT,NRECT, IOPTION,ICHECK, $ IBIG,NBIG, NZERO,NZLEFT, NZDOUBT1,NZDOUBT2, NZEROOLD, $ NELMAX,NPLINMAX,IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR, $ IASBACK,IDZHBACK,ICHKBACK, $ IPOSMODU,IPOSBOX, IPOS, NZEROQ, NZLEFTQ, JF, M, J, I, ISLASK, $ ISHEETLOK, ICUTTYP0OLD, ICWORKLIM0, $ IXROW, JCOLX, IEXPX, $ IAZ_TABNEW,IAZ_VALTERM,IAZ_TAB,IAZ_YTAB, MHQ(ntab), $ ISHEET(nzerodim), IZDOUBT(nzerodim), IARR(nzerodim), $ IARRLOK(nzerodim) REAL*8 FREQ,FREQMHQ0, HSTEPFSUBPRE, HSTEPFSUB0, XYTOL0,XYONENORM, $ FREQMIN,FREQINT, ALADAUER,PREALA,ALA, RIKTSORTANG, $ DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT, REDUCFAKT,DLOGREDUCFAKT,DDELT,XYTOL, $ UPILOK1,UPILOK2, UPIQQQ1,UPIQQQ2, RIKTSORTRE,RIKTSORTIM, $ YLOK,YLOKA,YLOKBCF, DLNDO,DLNDU, $ XYSLASK, SLASK, SLASK1,SLASK2, SLASKD, PI, $ OFREQ(nfrdim), UPREAL1(nrectdim),UPREAL2(nrectdim), $ UPIMAG1(nrectdim),UPIMAG2(nrectdim), ACOEF(jtab-1,ntab), $ HQI(jtab,2:ntab), FIHQI(jtab),QHQI(3), EPS0(1), XBIG(0:6), $ UPZRE(nzerodim), UPZIM(nzerodim), $ XYONEMIN(nzerodim), XYONEMAX(nzerodim), ARRLOK(nzerodim), $ UPZREQ(nzerodim), UPZIMQ(nzerodim), $ XYONEMINQ(nzerodim), XYONEMAXQ(nzerodim) COMPLEX*16 AZIMILN, UP,DERIVDO,DERIVDU, CSLASK, $ ZC(icworkdim) LOGICAL MORE3, OKSUB CHARACTER*32 FILM, MAINMODU,MAINBOX, FILMODU,FILBOX EXTERNAL CDETSUB,CDETSUBCIRC EXTERNAL IAZ_TABNEW,IAZ_VALTERM,IAZ_TAB,IAZ_YTAB DATA PI /3.141592653589793D0/ RW0 = (0.0D0,0.0D0) ! only mode search in rmodfaut IDET = 1 IDHTYP = 0 IHARKLIQ = 0 IRKSPEC = 0 WRITE(*,*) WRITE(*,*) 'Give infile FILM:' FILM = ' ' DO WHILE (FILM.EQ.' ') READ(*,'(A)') FILM END DO CALL BLANKTOEND (FILM) WRITE(*,*) 'Give MHQ0, FREQMHQ0 (Hz), IAZIMI :' READ(*,*) MHQ0,FREQMHQ0, IAZIMI IF (IAZIMI.NE.0) THEN WRITE(*,*) 'Give ALFA1AZIM:' READ(*,*) ALFA1AZIM END IF WRITE(*,*) 'Give ICUTTYP0 ( 1="axis-hyperbolic for up" 2=', $ '"vertical for up" ' WRITE(*,*) '3="mixture with aimcut" ) :' WRITE(*,*) '( "-5 or -15" for Muir )' WRITE(*,*) '( irrelevant for a rigid or free bottom )' READ(*,*) ICUTTYP0 ! used as a "work variable" IF (ICUTTYP0.EQ.3) THEN WRITE(*,*) 'Give AIMCUT:' READ(*,*) AIMCUT END IF IF (ICUTTYP0.NE.1 .AND. ICUTTYP0.NE.2 .AND. ICUTTYP0.NE.3 . $ AND. ICUTTYP0.NE.-5 .AND. ICUTTYP0.NE.-15) STOP 340 !Muirpoles notimpl ICUTTYP = ICUTTYP0 IF (ICUTTYP.EQ.1) ICUTTYP = 2 WRITE(*,*) 'Give ISHEETTYP0 ( 1="++" 2="-+" 3="+-" 4="--" ) for up" :' WRITE(*,*) '( irrelevant for a rigid or free bottom )' WRITE(*,*) '( "1 or 3" "2 or 4" irrelevant for a fluid half-space )' READ(*,*) ISHEETTYP0 WRITE(*,*) 'Give EPSGLOBLN,EPSGLOB,EPSFAKT,HFAKTMAX,HOFF,', $ 'EPSMAGNIFMAX :' WRITE(*,*) '( standard is epsfakt=0.1, hfaktmax=1.05, hoff=6d-16 )' READ(*,*) EPSGLOBLN,EPSGLOB,EPSFAKT,HFAKTMAX,HOFF,EPSMAGNIFMAX WRITE(*,*) 'Give EPSFSUB,HSTEPFSUBPRE :' WRITE(*,*) '( non-positive epsfsub implies "dpclin with ', $ 'constant hstep without ' WRITE(*,*) 'internal error control )' WRITE(*,*) '( except for fsmdzmin, hstepfsub will be ', $ '"hstepfsubpre/freq" which ' WRITE(*,*) ' is a "compromise" regarding rkra & rkrb )' READ(*,*) EPSFSUB,HSTEPFSUBPRE WRITE(*,*) 'Give (main) XYTOL0, IXTYP0,NRUS0,IRUS0,JRUS0 :' WRITE(*,*) '( xytol0 is an absolute tolerance )' WRITE(*,*) '( ixtyp0 = "1,-1,0 / 2,-2 / 3,-3" for ', $ '"polynomial / rational / inv-poly" )' READ(*,*) XYTOL0, IXTYP0,NRUS0,IRUS0,JRUS0 IF ( NRUS0.NE.0 .AND. IOPTION.EQ.0 ) THEN WRITE(*,*) ' %%% extrapolation not relevant' NRUS0 = 0 END IF IF (IABS(IXTYP0).GE.3) STOP 5710 IF ( NRUS0.EQ.1 .OR. NRUS0.GT.ntab ) STOP 571 IF (MIN0(NRUS0,JRUS0).GT.jtab) STOP 572 IF (NRUS0.GT.0) THEN WRITE(*,*) 'Give XYONENORM:' READ(*,*) XYONENORM END IF WRITE(*,*) WRITE(*,*) 'Give FREQMIN(>0.0),FREQINT (both in Hz) , NFR :' WRITE(*,*) '( " freqmin<0.0 & freqint=0.0 " for ', $ '"explicit ofreq(1:nfr) " )' READ(*,*) FREQMIN,FREQINT, NFR IF ( FREQMIN.EQ.0.0D0 .OR. NFR.GT.nfrdim ) STOP 6666 IF (FREQMIN.GE.0.0D0) THEN SLASK = FREQINT/NFR DO JF = 1, NFR OFREQ(JF) = FREQMIN + (JF-1)*SLASK END DO ELSE IF (FREQINT.NE.0.0D0) STOP 6660 WRITE(*,*) 'give ofreq(1:nfr) increasing' READ(*,*) (OFREQ(JF),JF=1,NFR) END IF WRITE(*,*) 'Give ALADAUER (s), PREALA :' WRITE(*,*) '( 0.0 ' DO IRECT = 1, NRECT READ(*,*) UPREAL1(IRECT),UPREAL2(IRECT),UPIMAG1(IRECT),UPIMAG2(IRECT) END DO WRITE(*,*) WRITE(*,*) 'Give IOPTION ( 0 for "just nzero", 1 for ', $ '"normal", 2 for "+ log" ) :' READ(*,*) IOPTION WRITE(*,*) 'Give ICHECK ( 1 for additional "argseg" check ) :' READ(*,*) ICHECK WRITE(*,*) 'Give RIKTSORTANG (in degrees) :' READ(*,*) RIKTSORTANG RIKTSORTRE = DCOSD(RIKTSORTANG) RIKTSORTIM = DSIND(RIKTSORTANG) WRITE(*,*) 'Give "main name" MAINMODU for the nfr "outfiles-mode" :' WRITE(*,*) '(freq indices as extensions, plot with plotfkn)' READ(*,'(A)') MAINMODU CALL BLANKTOEND (MAINMODU) I = 32 DO WHILE (MAINMODU(I:I).EQ.' ') I = I-1 END DO MAINMODU((I+1):(I+1)) = '.' IPOSMODU = I+2 WRITE(*,*) 'Give "main name" MAINBOX for the nfr ', $ '"outfiles-box(remaining)" :' WRITE(*,*) '(freq indices as extensions)' READ(*,'(A)') MAINBOX CALL BLANKTOEND (MAINBOX) I = 32 DO WHILE (MAINBOX(I:I).EQ.' ') I = I-1 END DO MAINBOX((I+1):(I+1)) = '.' IPOSBOX = I+2 WRITE(*,*) WRITE(*,*) 'Give DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX:' READ(*,*) DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX WRITE(*,*) 'Give DERIVDRELMAX,DRELMAX,DRELGOALFAKT:' READ(*,*) DERIVDRELMAX,DRELMAX,DRELGOALFAKT WRITE(*,*) 'Give ARGDIFSTORE:' READ(*,*) ARGDIFSTORE WRITE(*,*) 'Give SLIMLEFT:' READ(*,*) SLIMLEFT WRITE(*,*) 'Give REDUCFAKT,DDELT,XYTOL:' READ(*,*) REDUCFAKT,DDELT,XYTOL DLOGREDUCFAKT = DLOG(REDUCFAKT) OPEN(UNIT=20,FILE=FILM,STATUS='OLD') ISLASK = IAZIMI IAZIMI = 0 ! for this initial "particular" call CALL RMODEL ( 20,0, MHQ0 ) NSOURCFW = 0 ! irrelevant here IAZIMI = ISLASK ! reset iazimi CLOSE(20) IF ( NSW.EQ.0 .AND. (ISOLID(0).NE.0.OR.RHO(0).LT.0.0D0) ) IHARKLIQ = 1 C to be able to treat 'all-solid' in dhcomp (rigid or half-space top) IF (NRUS0.GT.0) THEN IF ( IRUS0.GT.NRUS0 .OR. JRUS0.LT.2 ) STOP 447 IADR00_TABNEW = 0 ! "tabnew(1:1)" <--> zc(1) IADR00_VALTERM = IADR00_TABNEW + 1 IADR00_TAB = IADR00_VALTERM + 1 ICWORKLIM0 = IADR00_TAB + MIN0(NRUS0,JRUS0)*1 IF (ICWORKLIM0.GT.icworkdim) STOP 451 C update the "Bulirsch Folge" MHQ(1) = 1 MHQ(2) = 2 MHQ(3) = 3 IXROW = 4 DO WHILE (IXROW.LE.NRUS0) MHQ(IXROW) = 2 * MHQ(IXROW-2) IXROW = IXROW + 2 END DO IXROW = 5 DO WHILE (IXROW.LE.NRUS0) MHQ(IXROW) = 2 * MHQ(IXROW-2) IXROW = IXROW + 2 END DO EPS0(1) = XYTOL0 C computation of "hqi(,)" IF ( IIFSMSOME.EQ.0 .OR. EPSFSUB.GT.0.0D0 ) THEN IF (IIMHQSOME.EQ.0) THEN WRITE(*,*) ' %%% extrapolation not needed' NRUS0 = 0 ELSE C now "I(h)" series with the even powers of h ("ODE trapezoidal") DO I=2,NRUS0 ! i=iminx,nrus JCOLX = MIN0(I,JRUS0) DO J=2,JCOLX HQI(J,I) = ( DBLE(MHQ(I)) / MHQ(I-J+1) ) ** 2 C ! "i"/"i-j+1" END DO END DO END IF ELSE IF (IABS(IXTYP0).NE.2) THEN JCOLX = MIN0(NRUS0,JRUS0) DO M=1,JCOLX-1 IF (IIMHQSOME.EQ.0) THEN C now "I(h)" series with h**5,h**6,... (Dormand-Prince global C error) IEXPX = 4+M ELSE C now "I(h)" series with h**2,h**4,h**5,h**6,... (combination) IF (M.EQ.1) THEN IEXPX = 2 ELSE IEXPX = 2+M END IF END IF DO I=1,NRUS0 ACOEF(M,I) = ( DBLE(MHQ(1)) / MHQ(I) ) ** IEXPX END DO END DO DO J=2,JCOLX DO I=NRUS0,J,-1 ! iminx=2 SLASK = ACOEF(J-1,I-1) / ACOEF(J-1,I) HQI(J,I) = SLASK SLASKD = 1 / (SLASK-1) DO M=J,JCOLX-1 ACOEF(M,I) = SLASKD * ( SLASK*ACOEF(M,I) - ACOEF(M,I-1) ) END DO END DO END DO ELSE C now "I(h)" series with the even powers of h, disregarding C a priori information on "gaps" in the asymptotic series DO I=2,NRUS0 ! i=iminx,nrus JCOLX = MIN0(I,JRUS0) DO J=2,JCOLX HQI(J,I) = ( DBLE(MHQ(I)) / MHQ(I-J+1) ) ** 2 C ! "i"/"i-j+1" END DO END DO END IF END IF DO JF = 1, NFR WRITE(*,*) WRITE(*,*) 'jf,nfr:',JF,NFR FREQ = OFREQ(JF) OMEG = DCMPLX ( 2*PI*FREQ, ALA ) IF (IAZIMI.NE.0) THEN CSLASK = ALFA1AZIM * OMEG FREAZIM = 2 * AZIMILN(CSLASK) / PI FIMAZIM = CSLASK / (1+CSLASK) END IF MHQNOW = FREQ/FREQMHQ0 MHQNOW = MHQ0 * JMAX0(1,MHQNOW) MHQEFF = MHQNOW CALL RMODEL ( 0,0, MHQEFF ) IF (IUROT.EQ.0) THEN U0 = (1.0D0,0.0D0) RK0 = OMEG ELSE RK0 = CDABS(OMEG) ! rk0 is now positive real U0 = RK0 / OMEG ! rk0=omeg*u0, arg(u0)=-arg(omeg) END IF OMEG2 = OMEG*OMEG IF (RHO(NOL).GT.0.0D0) THEN C note: the rare cases with pi/2 <= arg(upbpa) < pi or C pi/2 <= arg(upbpb) < pi , possible only for truly complex omeg, C were not considered when the "call zerofind" structure below was C programmed UPBPA = CDSQRT(UA2(NOL)) / U0 UPBPA2 = UPBPA * UPBPA IF (DREAL(UPBPA).LE.0.0D0) STOP 707 ! see above IF (ICUTTYP0.EQ.1) XYHYPA = DREAL(UPBPA) * DIMAG(UPBPA) WRITE(*,*) WRITE(*,*) 'upbpa is ',UPBPA IF (ISOLID(NOL).NE.0) THEN UPBPB = CDSQRT(UB2(NOL)) / U0 UPBPB2 = UPBPB * UPBPB IF (DREAL(UPBPB).LE.DREAL(UPBPA)) STOP 707 ! rare case not impl IF (ICUTTYP0.EQ.1) XYHYPB = DREAL(UPBPB) * DIMAG(UPBPB) WRITE(*,*) 'upbpb is ',UPBPB ELSE UPBPB = (0.0D0,0.0D0) UPBPB2 = (0.0D0,0.0D0) END IF ELSE UPBPA = (0.0D0,0.0D0) UPBPA2 = (0.0D0,0.0D0) UPBPB = (0.0D0,0.0D0) UPBPB2 = (0.0D0,0.0D0) END IF IF (RHO(0).GT.0.0D0) THEN C note: the rare cases with pi/2 <= arg(upbpasurf) < pi or C pi/2 <= arg(upbpbsurf) < pi , possible only for truly cmplx omeg, C were not considered when the "call zerofind" structure below was C programmed UPBPASURF = CDSQRT(UA2(0)) / U0 UPBPA2SURF = UPBPASURF * UPBPASURF IF (DREAL(UPBPASURF).LE.0.0D0) STOP 707 ! see above IF (ICUTTYP0.NE.2) THEN WRITE(*,*) 'upper half-space, vertical branch-cuts needed' STOP 341552 END IF WRITE(*,*) WRITE(*,*) 'upbpasurf is ',UPBPASURF IF (ISOLID(0).NE.0) THEN UPBPBSURF = CDSQRT(UB2(0)) / U0 UPBPB2SURF = UPBPBSURF * UPBPBSURF IF (DREAL(UPBPBSURF).LE.DREAL(UPBPASURF)) STOP 707 !rare, no impl IF (ICUTTYP0.NE.2) STOP 341552 WRITE(*,*) 'upbpbsurf is ',UPBPBSURF ELSE UPBPBSURF = (0.0D0,0.0D0) UPBPB2SURF = (0.0D0,0.0D0) END IF IF (ICUTTYP0.NE.2) STOP 341552 ELSE UPBPASURF = (0.0D0,0.0D0) UPBPA2SURF = (0.0D0,0.0D0) UPBPBSURF = (0.0D0,0.0D0) UPBPB2SURF = (0.0D0,0.0D0) END IF FILMODU = MAINMODU IPOS = IPOSMODU ISLASK = FREQ CALL PUTINTP (ISLASK,IPOS,FILMODU) IF (ISLASK.NE.FREQ) THEN IF (IPOS.GT.32) STOP 460 FILMODU(IPOS:IPOS) = 'X' IPOS = IPOS+1 CALL PUTINTP (JF,IPOS,FILMODU) END IF OPEN(UNIT=1,FILE=FILMODU,STATUS='NEW') ISLASK = 24 + 3*NRECT IF (IAZIMI.NE.0) ISLASK = ISLASK+1 IF (ICUTTYP0.EQ.3) ISLASK = ISLASK+1 IF (NRUS0.GT.0) ISLASK = ISLASK+1 WRITE(1,*) ISLASK,SNGL(FREQ),' ! "lines to skip", freq' WRITE(1,*) ' FILMODU from RMODFAUT ; freq=',SNGL(FREQ) WRITE(1,*) ' FILM: ',FILM WRITE(1,*) ' MHQ0,FREQMHQ0, IAZIMI: ',MHQ0,SNGL(FREQMHQ0), IAZIMI IF (IAZIMI.NE.0) WRITE(1,*) ' ALFA1AZIM: ',SNGL(ALFA1AZIM) IF (ICUTTYP0.NE.3) THEN WRITE(1,*) ' ICUTTYP0,ISHEETTYP0: ',ICUTTYP0,ISHEETTYP0 ELSE WRITE(1,*) ' ICUTTYP0,AIMCUT: ',ICUTTYP0,SNGL(AIMCUT) WRITE(1,*) ' ISHEETTYP0: ',ISHEETTYP0 END IF WRITE(1,*) ' EPSGLOBLN,EPSGLOB: ',SNGL(EPSGLOBLN),SNGL(EPSGLOB) WRITE(1,*) ' EPSFAKT,HFAKTMAX: ',SNGL(EPSFAKT),SNGL(HFAKTMAX) WRITE(1,*) ' HOFF,EPSMAGNIFMAX: ',SNGL(HOFF),SNGL(EPSMAGNIFMAX) WRITE(1,*) ' EPSFSUB,HSTEPFSUBPRE: ',SNGL(EPSFSUB), $ SNGL(HSTEPFSUBPRE) WRITE(1,*) ' XYTOL0, IXTYP0 : ',SNGL(XYTOL0),IXTYP0 WRITE(1,*) ' NRUS0,IRUS0,JRUS0 : ',NRUS0,IRUS0,JRUS0 IF (NRUS0.GT.0) WRITE(1,*) ' XYONENORM : ',XYONENORM WRITE(1,*) ' FREQ: ',FREQ WRITE(1,*) ' ALADAUER,PREALA: ',SNGL(ALADAUER),SNGL(PREALA) WRITE(1,*) ' IUROT: ',IUROT WRITE(1,*) ' NRECT: ',NRECT DO IRECT = 1, NRECT WRITE(1,*) ' for irect=',IRECT WRITE(1,*) ' UPREAL1&2 ',SNGL(UPREAL1(IRECT)),SNGL(UPREAL2(IRECT)) WRITE(1,*) ' UPIMAG1&2 ',SNGL(UPIMAG1(IRECT)),SNGL(UPIMAG2(IRECT)) END DO WRITE(1,*) ' IOPTION,ICHECK: ',IOPTION,ICHECK WRITE(1,*) ' RIKTSORTANG: ',RIKTSORTANG WRITE(1,*) ' MAINMODU ',MAINMODU WRITE(1,*) ' MAINBOX ',MAINBOX WRITE(1,*) ' DDSUPMIN,DDSUPMAX: ',SNGL(DDSUPMIN),SNGL(DDSUPMAX) WRITE(1,*) ' DDSUPFAKTMAX: ',SNGL(DDSUPFAKTMAX) WRITE(1,*) ' DERIVDRELMAX,DRELMAX: ',SNGL(DERIVDRELMAX),SNGL(DRELMAX) WRITE(1,*) ' DRELGOALFAKT, ARGDIFSTORE: ', $ SNGL(DRELGOALFAKT),' ',SNGL(ARGDIFSTORE) WRITE(1,*) ' SLIMLEFT: ',SNGL(SLIMLEFT) WRITE(1,*) ' REDUCFAKT,DDELT,XYTOL: ',SNGL(REDUCFAKT), $ SNGL(DDELT),SNGL(XYTOL) WRITE(1,*) IF (NRUS0.LE.0) THEN WRITE(1,*) 4, ICUTTYP0,ISHEETTYP0, ALA,IUROT ELSE WRITE(1,*) 5, ICUTTYP0,ISHEETTYP0, ALA,IUROT END IF WRITE(1,*) FILBOX = MAINBOX IPOS = IPOSBOX ISLASK = FREQ CALL PUTINTP (ISLASK,IPOS,FILBOX) IF (ISLASK.NE.FREQ) THEN IF (IPOS.GT.32) STOP 460 FILBOX(IPOS:IPOS) = 'X' IPOS = IPOS+1 CALL PUTINTP (JF,IPOS,FILBOX) END IF OPEN(UNIT=11,FILE=FILBOX,STATUS='NEW') WRITE(11,*) ' FILBOX from RMODFAUT ; freq=',SNGL(FREQ) WRITE(11,*) ' MAINMODU ',MAINMODU WRITE(11,*) ' MAINBOX ',MAINBOX WRITE(11,*) IF (IOPTION.EQ.2) THEN OPEN(UNIT=99,STATUS='NEW') WRITE(99,*) ' "zerofind log" from rmodfaut ; freq=',SNGL(FREQ) WRITE(99,*) ' MAINMODU ',MAINMODU WRITE(99,*) ' MAINBOX ',MAINBOX END IF IF (IIFSMSOME.NE.0) THEN IF (EPSFSUB.GT.0.0D0) THEN HSTEPFSUB = HSTEPFSUBPRE/FREQ ELSE HSTEPFSUB0 = DMIN1 ( HSTEPFSUBPRE/FREQ, FSMDZMIN ) ! to avoid C extrapolating on "equal hstepfsub" HSTEPFSUB = HSTEPFSUB0 END IF END IF C ------------------------------------------------------------------ C main loop (beginning after some initialization) NZERO = 0 NZLEFT = 0 NELMAX = 0 NPLINMAX = 0 IASFCALL = 0 IDZHFCALL = 0 ICHKFCALL = 0 INELCOMPR = 0 IASBACK = 0 IDZHBACK = 0 ICHKBACK = 0 DO IRECT = 1, NRECT UPILOK1 = UPIMAG1(IRECT) UPILOK2 = UPIMAG2(IRECT) C note: C *** with the "nbig split" for vertical branch-cuts, parts of the C "inner" vertical boundaries will be treated twice; this could C of course be avoided by basing the "nbig split" on "y values" as C well (and not only "x values") C *** with repeated computations for different "isheettyp", as may be C done below when icuttyp0=1, the "inner" branch-cuts are treated C more than once which could also be avoided in principle, but this C is probably not worth while C *** the "call zerofind" structure below for "icuttyp0=1 & upbpa C nonzero" is by no means optimal NBIG = 1 XBIG(0) = UPREAL1(IRECT) XBIG(1) = UPREAL2(IRECT) IF ( UPBPA.NE.(0.0D0,0.0D0) .OR. UPBPASURF.NE.(0.0D0,0.0D0) ) THEN IF (UPBPA.NE.(0.0D0,0.0D0)) THEN SLASK = DIMAG(UPBPA) ELSE IF (UPBPASURF.NE.(0.0D0,0.0D0)) THEN SLASK = DIMAG(UPBPASURF) END IF IF (UPBPA.NE.(0.0D0,0.0D0)) SLASK = DMIN1(SLASK,DIMAG(UPBPA)) IF (UPBPB.NE.(0.0D0,0.0D0)) SLASK = DMIN1(SLASK,DIMAG(UPBPB)) IF (UPBPASURF.NE.(0.0D0,0.0D0)) SLASK = $ DMIN1(SLASK,DIMAG(UPBPASURF)) IF (UPBPBSURF.NE.(0.0D0,0.0D0)) SLASK = $ DMIN1(SLASK,DIMAG(UPBPBSURF)) IF ( ( XBIG(1).LE.0.0D0 .AND. UPILOK1.GT.-SLASK ) .OR. $ ( XBIG(0).GE.0.0D0 .AND. UPILOK2.LT.SLASK ) ) THEN IF ( ICUTTYP0.EQ.1 .OR. ICUTTYP0.EQ.3 ) ICUTTYP0 = 2 ! possible GOTO 141 ! simple case END IF END IF DO I = 1, 9 IF (I.LE.2) THEN IF (UPBPA.NE.(0.0D0,0.0D0)) THEN IF (I.EQ.1) THEN SLASK = -DREAL(UPBPA) ELSE SLASK = DREAL(UPBPA) END IF ELSE GOTO 41 END IF ELSE IF (I.LE.4) THEN IF (UPBPB.NE.(0.0D0,0.0D0)) THEN IF (I.EQ.3) THEN SLASK = -DREAL(UPBPB) ELSE SLASK = DREAL(UPBPB) END IF ELSE GOTO 41 END IF ELSE IF (I.EQ.5) THEN IF ( (ICUTTYP0.EQ.1.OR.ICUTTYP0.EQ.3) .AND. $ ( UPBPA.NE.(0.0D0,0.0D0).OR.UPBPASURF.NE.(0.0D0,0.0D0) ) ) THEN SLASK = 0.0D0 ELSE GOTO 41 END IF ELSE IF (I.LE.7) THEN IF (UPBPASURF.NE.(0.0D0,0.0D0)) THEN IF (I.EQ.6) THEN SLASK = -DREAL(UPBPASURF) ELSE SLASK = DREAL(UPBPASURF) END IF ELSE GOTO 41 END IF ELSE IF (UPBPBSURF.NE.(0.0D0,0.0D0)) THEN IF (I.EQ.8) THEN SLASK = -DREAL(UPBPBSURF) ELSE SLASK = DREAL(UPBPBSURF) END IF ELSE GOTO 41 END IF END IF IF (SLASK.GT.XBIG(0)) THEN IF (SLASK.LT.XBIG(NBIG)) THEN ISLASK = 1 DO WHILE (SLASK.GT.XBIG(ISLASK)) ISLASK = ISLASK + 1 END DO IF (XBIG(ISLASK).GT.SLASK) THEN DO IBIG = NBIG, ISLASK, -1 XBIG(IBIG+1) = XBIG(IBIG) END DO XBIG(ISLASK) = SLASK NBIG = NBIG + 1 END IF END IF END IF 41 CONTINUE END DO 141 WRITE(*,*) WRITE(*,*) '"irect,nbig" now ',IRECT,NBIG WRITE(*,*) DO IBIG = 1, NBIG UPREF = DCMPLX ( 0.5D0*(XBIG(IBIG-1)+XBIG(IBIG)) , $ 0.5D0*(UPILOK1+UPILOK2) ) ICUTTYP0OLD = ICUTTYP0 IF (ICUTTYP0.EQ.3) THEN C sometimes the "split" that may follow is unnecessary IF (DREAL(UPREF).GT.0.0D0) THEN IF (UPILOK2.LE.AIMCUT) THEN ICUTTYP0 = 1 MORE3 = .FALSE. ELSE IF (AIMCUT.LE.UPILOK1) THEN ICUTTYP0 = 2 MORE3 = .FALSE. ELSE ICUTTYP0 = 1 UPILOK2 = AIMCUT MORE3 = .TRUE. END IF ELSE IF (UPILOK2.LE.-AIMCUT) THEN ICUTTYP0 = 2 MORE3 = .FALSE. ELSE IF (-AIMCUT.LE.UPILOK1) THEN ICUTTYP0 = 1 MORE3 = .FALSE. ELSE ICUTTYP0 = 2 UPILOK2 = -AIMCUT MORE3 = .TRUE. END IF END IF END IF GOTO 71 70 IF (ICUTTYP0.EQ.1) THEN ICUTTYP0 = 2 ELSE ICUTTYP0 = 1 END IF UPILOK1 = UPILOK2 UPILOK2 = UPIMAG2(IRECT) MORE3 = .FALSE. 71 CONTINUE IF ( ICUTTYP0.NE.1 .OR. UPBPA.EQ.(0.0D0,0.0D0) ) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROOLD = NZERO ICORNTEST = 0 CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPILOK1,UPILOK2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, ARGDIFSTORE, $ SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, IOPTION,ICHECK,11, $ nzerodim, NZERO,NZLEFT, UPZRE,UPZIM, XYONEMIN,XYONEMAX, $ NELMAX,NPLINMAX,IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR, $ IASBACK,IDZHBACK,ICHKBACK) DO I = NZEROOLD+1, NZERO ISHEET(I) = ISHEETTYP IZDOUBT(I) = 0 END DO ELSE IF (UPBPB.EQ.(0.0D0,0.0D0)) THEN IF ( DABS(DREAL(UPREF)) .GT. DREAL(UPBPA) ) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROOLD = NZERO ICORNTEST = 0 CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPILOK1,UPILOK2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, ARGDIFSTORE, $ SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, IOPTION,ICHECK,11, $ nzerodim, NZERO,NZLEFT, UPZRE,UPZIM, XYONEMIN,XYONEMAX, $ NELMAX,NPLINMAX,IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR, $ IASBACK,IDZHBACK,ICHKBACK) DO I = NZEROOLD+1, NZERO ISHEET(I) = ISHEETTYP IZDOUBT(I) = 0 END DO ELSE IF (DREAL(UPREF).GT.0.0D0) THEN IF (XBIG(IBIG).EQ.DREAL(UPBPA)) THEN YLOK = DIMAG(UPBPA) - DDSUPMAX !avoid "argseg" start at upbpa ELSE YLOK = XYHYPA / XBIG(IBIG) END IF IF (YLOK.EQ.0.0D0) YLOK = -DDSUPMAX !avoid the axis IF (UPILOK1.LT.YLOK) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROOLD = NZERO ICORNTEST = 0 UPIQQQ2 = DMIN1(YLOK,UPILOK2) CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPILOK1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZERO,NZLEFT, UPZRE,UPZIM, $ XYONEMIN,XYONEMAX, NELMAX,NPLINMAX,IASFCALL,IDZHFCALL, $ ICHKFCALL,INELCOMPR,IASBACK,IDZHBACK,ICHKBACK) DO I = NZEROOLD+1, NZERO ISHEET(I) = ISHEETTYP IZDOUBT(I) = 0 END DO END IF IF (YLOK.LT.UPILOK2) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 1 UPIQQQ1 = DMAX1(YLOK,UPILOK1) IF (XBIG(IBIG-1).GT.0.0D0) THEN UPIQQQ2 = XYHYPA/XBIG(IBIG-1) IF (UPIQQQ2.EQ.0.0D0) UPIQQQ2 = +DDSUPMAX !avoid the axis UPIQQQ2 = DMIN1 ( UPILOK2, UPIQQQ2 ) ELSE UPIQQQ2 = UPILOK2 END IF IF (UPIQQQ1.LT.UPIQQQ2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, UPZREQ, $ UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .LE. XYHYPA ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. XYHYPA ) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. XYHYPA ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE.0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. XYHYPA ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. XYHYPA ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF ISHEETTYP = ISHEETLOK(1,ISHEETTYP0) WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 2 UPIQQQ1 = DMAX1 ( YLOK, UPILOK1 ) IF (UPIQQQ1.LT.UPILOK2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPILOK2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, UPZREQ, $ UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .GE. XYHYPA ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE.0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. XYHYPA ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. XYHYPA ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF END IF ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. XYHYPA ) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. XYHYPA ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF END IF ELSE IF (XBIG(IBIG-1).EQ.DREAL(-UPBPA)) THEN YLOK = DIMAG(-UPBPA) + DDSUPMAX !avoi "argseg" strt at -upbpa ELSE YLOK = XYHYPA / XBIG(IBIG-1) END IF IF (YLOK.EQ.0.0D0) YLOK = +DDSUPMAX !avoid the axis IF (YLOK.LT.UPILOK2) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROOLD = NZERO ICORNTEST = 0 UPIQQQ1 = DMAX1(YLOK,UPILOK1) CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPILOK2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZERO,NZLEFT, UPZRE, $ UPZIM, XYONEMIN,XYONEMAX, NELMAX,NPLINMAX,IASFCALL, $ IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK,IDZHBACK, $ ICHKBACK) DO I = NZEROOLD+1, NZERO ISHEET(I) = ISHEETTYP IZDOUBT(I) = 0 END DO END IF IF (UPILOK1.LT.YLOK) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 1 IF (XBIG(IBIG).LT.0.0D0) THEN UPIQQQ1 = XYHYPA/XBIG(IBIG) IF (UPIQQQ1.EQ.0.0D0) UPIQQQ1 = -DDSUPMAX !avoid the axis UPIQQQ1 = DMAX1 ( UPILOK1, UPIQQQ1 ) ELSE UPIQQQ1 = UPILOK1 END IF UPIQQQ2 = DMIN1(YLOK,UPILOK2) IF (UPIQQQ1.LT.UPIQQQ2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, UPZREQ, $ UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .LE. XYHYPA ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. XYHYPA ) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. XYHYPA ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE.0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. XYHYPA ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. XYHYPA ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF ISHEETTYP = ISHEETLOK(1,ISHEETTYP0) WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 2 UPIQQQ2 = DMIN1 ( YLOK, UPILOK2 ) IF (UPILOK1.LT.UPIQQQ2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPILOK1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, UPZREQ, $ UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .GE. XYHYPA ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE.0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. XYHYPA ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. XYHYPA ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF END IF ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. XYHYPA ) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. XYHYPA ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF END IF END IF ELSE IF ( DABS(DREAL(UPREF)) .GT. DREAL(UPBPB) ) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROOLD = NZERO ICORNTEST = 0 CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPILOK1,UPILOK2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, ARGDIFSTORE, $ SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, IOPTION,ICHECK,11, $ nzerodim, NZERO,NZLEFT, UPZRE,UPZIM, XYONEMIN, $ XYONEMAX, NELMAX,NPLINMAX,IASFCALL,IDZHFCALL, $ ICHKFCALL,INELCOMPR,IASBACK,IDZHBACK,ICHKBACK) DO I = NZEROOLD+1, NZERO ISHEET(I) = ISHEETTYP IZDOUBT(I) = 0 END DO ELSE IF ( DABS(DREAL(UPREF)) .GT. DREAL(UPBPA) ) THEN IF (DREAL(UPREF).GT.0.0D0) THEN IF (XBIG(IBIG).EQ.DREAL(UPBPB)) THEN YLOK = DIMAG(UPBPB) - DDSUPMAX !avoi "argseg" strt at upbpb ELSE YLOK = XYHYPB / XBIG(IBIG) END IF IF (YLOK.EQ.0.0D0) YLOK = -DDSUPMAX !avoid the axis IF (UPILOK1.LT.YLOK) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROOLD = NZERO ICORNTEST = 0 UPIQQQ2 = DMIN1(YLOK,UPILOK2) CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPILOK1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZERO,NZLEFT, UPZRE, $ UPZIM, XYONEMIN,XYONEMAX, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = NZEROOLD+1, NZERO ISHEET(I) = ISHEETTYP IZDOUBT(I) = 0 END DO END IF IF (YLOK.LT.UPILOK2) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 3 UPIQQQ1 = DMAX1(YLOK,UPILOK1) IF (XBIG(IBIG-1).GT.0.0D0) THEN UPIQQQ2 = XYHYPB/XBIG(IBIG-1) IF (UPIQQQ2.EQ.0.0D0) UPIQQQ2 = +DDSUPMAX !avoid the axis UPIQQQ2 = DMIN1 ( UPILOK2, UPIQQQ2 ) ELSE UPIQQQ2 = UPILOK2 END IF IF (UPIQQQ1.LT.UPIQQQ2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, UPZREQ, $ UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .LE. XYHYPB ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. XYHYPB ) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. XYHYPB ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE. $ 0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. XYHYPB ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. XYHYPB ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF ISHEETTYP = ISHEETLOK(2,ISHEETTYP0) WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 4 UPIQQQ1 = DMAX1 ( YLOK, UPILOK1 ) IF (UPIQQQ1.LT.UPILOK2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPILOK2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, UPZREQ, $ UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .GE. XYHYPB ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE. $ 0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. XYHYPB ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. XYHYPB ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF END IF ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. XYHYPB ) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. XYHYPB ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF END IF ELSE IF (XBIG(IBIG-1).EQ.DREAL(-UPBPB)) THEN YLOK = DIMAG(-UPBPB) + DDSUPMAX !avoi"argseg"strt at -upbpb ELSE YLOK = XYHYPB / XBIG(IBIG-1) END IF IF (YLOK.EQ.0.0D0) YLOK = +DDSUPMAX !avoid the axis IF (YLOK.LT.UPILOK2) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROOLD = NZERO ICORNTEST = 0 UPIQQQ1 = DMAX1(YLOK,UPILOK1) CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPILOK2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZERO,NZLEFT, UPZRE, $ UPZIM, XYONEMIN,XYONEMAX, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = NZEROOLD+1, NZERO ISHEET(I) = ISHEETTYP IZDOUBT(I) = 0 END DO END IF IF (UPILOK1.LT.YLOK) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 3 IF (XBIG(IBIG).LT.0.0D0) THEN UPIQQQ1 = XYHYPB/XBIG(IBIG) IF (UPIQQQ1.EQ.0.0D0) UPIQQQ1 = -DDSUPMAX !avoid the axis UPIQQQ1 = DMAX1 ( UPILOK1, UPIQQQ1 ) ELSE UPIQQQ1 = UPILOK1 END IF UPIQQQ2 = DMIN1(YLOK,UPILOK2) IF (UPIQQQ1.LT.UPIQQQ2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, UPZREQ, $ UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .LE. XYHYPB ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. XYHYPB ) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. XYHYPB ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE. $ 0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. XYHYPB ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. XYHYPB ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF ISHEETTYP = ISHEETLOK(2,ISHEETTYP0) WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 4 UPIQQQ2 = DMIN1 ( YLOK, UPILOK2 ) IF (UPILOK1.LT.UPIQQQ2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPILOK1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, UPZREQ, $ UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .GE. XYHYPB ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE. $ 0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. XYHYPB ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. XYHYPB ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF END IF ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. XYHYPB ) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. XYHYPB ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF END IF END IF ELSE IF (DREAL(UPREF).GT.0.0D0) THEN IF (XBIG(IBIG).EQ.DREAL(UPBPA)) THEN YLOKA = DIMAG(UPBPA) - DDSUPMAX !avoi"argseg"strt at upbpa ELSE YLOKA = XYHYPA / XBIG(IBIG) END IF IF (YLOKA.EQ.0.0D0) YLOKA = -DDSUPMAX !avoid the axis YLOKBCF = XYHYPB / XBIG(IBIG) IF (YLOKBCF.EQ.0.0D0) YLOKBCF = -DDSUPMAX !avoid the axis YLOK = DMIN1(YLOKA,YLOKBCF) IF (UPILOK1.LT.YLOK) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROOLD = NZERO ICORNTEST = 0 UPIQQQ2 = DMIN1(YLOK,UPILOK2) CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPILOK1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZERO,NZLEFT, UPZRE, $ UPZIM, XYONEMIN,XYONEMAX, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR, $ IASBACK,IDZHBACK,ICHKBACK) DO I = NZEROOLD+1, NZERO ISHEET(I) = ISHEETTYP IZDOUBT(I) = 0 END DO END IF IF (YLOK.LT.UPILOK2) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 5 SLASK = DMIN1(XYHYPA,XYHYPB) UPIQQQ1 = DMAX1(YLOK,UPILOK1) IF (XBIG(IBIG-1).GT.0.0D0) THEN UPIQQQ2 = SLASK/XBIG(IBIG-1) IF (UPIQQQ2.EQ.0.0D0) UPIQQQ2 = +DDSUPMAX !avoid the axis UPIQQQ2 = DMIN1 ( UPILOK2, UPIQQQ2 ) ELSE UPIQQQ2 = UPILOK2 END IF IF (UPIQQQ1.LT.UPIQQQ2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, UPZREQ, $ UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .LE. SLASK ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. SLASK ) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. SLASK ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE. $ 0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. SLASK ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. SLASK ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF ISHEETTYP = ISHEETLOK(3,ISHEETTYP0) WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 6 SLASK = DMAX1(XYHYPA,XYHYPB) UPIQQQ1 = DMAX1 ( YLOKA, YLOKBCF, UPILOK1 ) IF (UPIQQQ1.LT.UPILOK2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPILOK2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, $ UPZREQ,UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX, $ NPLINMAX,IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR, $ IASBACK,IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .GE. SLASK ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE. $ 0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. SLASK ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. SLASK ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF END IF ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. SLASK ) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. SLASK ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF IF (YLOKA.NE.YLOKBCF) THEN IF (YLOKA.LT.YLOKBCF) THEN ISHEETTYP = ISHEETLOK(1,ISHEETTYP0) ELSE ISHEETTYP = ISHEETLOK(2,ISHEETTYP0) END IF WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 7 SLASK1 = DMIN1(XYHYPA,XYHYPB) SLASK2 = DMAX1(XYHYPA,XYHYPB) UPIQQQ1 = DMAX1 ( YLOK, UPILOK1 ) IF (XBIG(IBIG-1).GT.0.0D0) THEN UPIQQQ2 = SLASK2/XBIG(IBIG-1) IF (UPIQQQ2.EQ.0.0D0) UPIQQQ2 = +DDSUPMAX !avoid axis UPIQQQ2 = DMIN1 ( UPILOK2, UPIQQQ2 ) ELSE UPIQQQ2 = UPILOK2 END IF IF (UPIQQQ1.LT.UPIQQQ2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, $ UPZREQ,UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX, $ NPLINMAX,IASFCALL,IDZHFCALL,ICHKFCALL, $ INELCOMPR,IASBACK,IDZHBACK,ICHKBACK) DO I = 1, NZEROQ SLASK = UPZREQ(I)*UPZIMQ(I) IF ( SLASK1.LE.SLASK .AND. SLASK.LE.SLASK2 ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. SLASK2 ) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. SLASK2 ) THEN IZDOUBT(NZERO) = 1 ELSE IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I). $ GE.0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. SLASK1 ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. SLASK1 ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF END IF END IF ELSE IF (SLASK2.LT.SLASK) THEN IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE. $ 0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. SLASK2 ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. SLASK2 ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. SLASK1 ) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. SLASK1 ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF END IF END IF ELSE IF (XBIG(IBIG-1).EQ.DREAL(-UPBPA)) THEN YLOKA = DIMAG(-UPBPA) +DDSUPMAX !avoi"argseg"strt at -upbpa ELSE YLOKA = XYHYPA / XBIG(IBIG-1) END IF IF (YLOKA.EQ.0.0D0) YLOKA = +DDSUPMAX !avoid the axis YLOKBCF = XYHYPB / XBIG(IBIG-1) IF (YLOKBCF.EQ.0.0D0) YLOKBCF = +DDSUPMAX !avoid the axis YLOK = DMAX1(YLOKA,YLOKBCF) IF (YLOK.LT.UPILOK2) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROOLD = NZERO ICORNTEST = 0 UPIQQQ1 = DMAX1(YLOK,UPILOK1) CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPILOK2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZERO,NZLEFT, UPZRE, $ UPZIM, XYONEMIN,XYONEMAX, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR, $ IASBACK,IDZHBACK,ICHKBACK) DO I = NZEROOLD+1, NZERO ISHEET(I) = ISHEETTYP IZDOUBT(I) = 0 END DO END IF IF (UPILOK1.LT.YLOK) THEN ISHEETTYP = ISHEETTYP0 WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 5 SLASK = DMIN1(XYHYPA,XYHYPB) IF (XBIG(IBIG).LT.0.0D0) THEN UPIQQQ1 = SLASK/XBIG(IBIG) IF (UPIQQQ1.EQ.0.0D0) UPIQQQ1 = -DDSUPMAX !avoid the axis UPIQQQ1 = DMAX1 ( UPILOK1, UPIQQQ1 ) ELSE UPIQQQ1 = UPILOK1 END IF UPIQQQ2 = DMIN1(YLOK,UPILOK2) IF (UPIQQQ1.LT.UPIQQQ2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, UPZREQ, $ UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR, $ IASBACK,IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .LE. SLASK ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. SLASK ) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. SLASK ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE. $ 0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. SLASK ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. SLASK ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF ISHEETTYP = ISHEETLOK(3,ISHEETTYP0) WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 6 SLASK = DMAX1(XYHYPA,XYHYPB) UPIQQQ2 = DMIN1 ( YLOKA, YLOKBCF, UPILOK2 ) IF (UPILOK1.LT.UPIQQQ2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPILOK1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, UPZREQ, $ UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX, $ IASFCALL,IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK, $ IDZHBACK,ICHKBACK) DO I = 1, NZEROQ IF ( (UPZREQ(I)*UPZIMQ(I)) .GE. SLASK ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE. $ 0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. SLASK ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. SLASK ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF END IF ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. SLASK ) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. SLASK ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF IF (YLOKA.NE.YLOKBCF) THEN IF (YLOKA.GT.YLOKBCF) THEN ISHEETTYP = ISHEETLOK(1,ISHEETTYP0) ELSE ISHEETTYP = ISHEETLOK(2,ISHEETTYP0) END IF WRITE(11,*) WRITE(11,*) ' --- ISHEETTYP NOW ',ISHEETTYP WRITE(11,*) NZEROQ = 0 ICORNTEST = 7 SLASK1 = DMIN1(XYHYPA,XYHYPB) SLASK2 = DMAX1(XYHYPA,XYHYPB) IF (XBIG(IBIG).LT.0.0D0) THEN UPIQQQ1 = SLASK2/XBIG(IBIG) IF (UPIQQQ1.EQ.0.0D0) UPIQQQ1 = -DDSUPMAX !avoid axis UPIQQQ1 = DMAX1 ( UPILOK1, UPIQQQ1 ) ELSE UPIQQQ1 = UPILOK1 END IF UPIQQQ2 = DMIN1 ( YLOK, UPILOK2 ) IF (UPIQQQ1.LT.UPIQQQ2) THEN CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ XBIG(IBIG-1),XBIG(IBIG), $ UPIQQQ1,UPIQQQ2, DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, $ DERIVDRELMAX,DRELMAX,DRELGOALFAKT, $ ARGDIFSTORE, SLIMLEFT,DLOGREDUCFAKT,DDELT,XYTOL, $ IOPTION,ICHECK,11,nzerodim, NZEROQ,NZLEFT, $ UPZREQ,UPZIMQ, XYONEMINQ,XYONEMAXQ, NELMAX, $ NPLINMAX,IASFCALL,IDZHFCALL,ICHKFCALL, $ INELCOMPR,IASBACK,IDZHBACK,ICHKBACK) DO I = 1, NZEROQ SLASK = UPZREQ(I)*UPZIMQ(I) IF ( SLASK1.LE.SLASK .AND. SLASK.LE.SLASK2 ) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. SLASK2 ) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. SLASK2 ) THEN IZDOUBT(NZERO) = 1 ELSE IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE. $ 0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. SLASK1 ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN IZDOUBT(NZERO) = 1 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. SLASK1 ) THEN IZDOUBT(NZERO) = 1 ELSE IZDOUBT(NZERO) = 0 END IF END IF END IF ELSE IF (SLASK2.LT.SLASK) THEN IF ( UPZREQ(I).GE.0.0D0 .AND. UPZIMQ(I).GE. $ 0.0D0 ) THEN IF (UPZREQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).LT.XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .LT. SLASK2 ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF ELSE IF (UPZREQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF (UPZIMQ(I).GT.-XYONEMINQ(I)) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .LT. SLASK2 ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF ELSE IF ( ((UPZREQ(I)+XYONEMINQ(I))* $ (UPZIMQ(I)+XYONEMINQ(I))) .GT. SLASK1 ) THEN ISLASK = 2 ELSE IF ( ((UPZREQ(I)-XYONEMINQ(I))* $ (UPZIMQ(I)-XYONEMINQ(I))) .GT. SLASK1 ) THEN ISLASK = 2 ELSE ISLASK = 0 END IF IF (ISLASK.NE.0) THEN IF (NZERO.GE.nzerodim) STOP 406 NZERO = NZERO + 1 UPZRE(NZERO) = UPZREQ(I) UPZIM(NZERO) = UPZIMQ(I) XYONEMIN(NZERO) = XYONEMINQ(I) XYONEMAX(NZERO) = XYONEMAXQ(I) ISHEET(NZERO) = ISHEETTYP IZDOUBT(NZERO) = 2 END IF END IF END DO END IF END IF END IF END IF END IF END IF IF (ICUTTYP0OLD.EQ.3) THEN IF (MORE3) THEN MORE3 = .FALSE. GOTO 70 ELSE ICUTTYP0 = ICUTTYP0OLD END IF END IF END DO END DO C end main loop C ------------------------------------------------------------------ C write results, sorted zeros in particular NZDOUBT1 = 0 NZDOUBT2 = 0 DO I = 1, NZERO IF (IZDOUBT(I).EQ.1) NZDOUBT1 = NZDOUBT1 + 1 IF (IZDOUBT(I).EQ.2) NZDOUBT2 = NZDOUBT2 + 1 END DO DO I = 1, NZERO ARRLOK(I) = UPZRE(I)*RIKTSORTRE + UPZIM(I)*RIKTSORTIM END DO CALL QUICKXD (ARRLOK,IARR,NZERO) DO I = 1, NZERO ARRLOK(I) = UPZRE(I) END DO DO I = 1, NZERO UPZRE(I) = ARRLOK(IARR(I)) END DO DO I = 1, NZERO ARRLOK(I) = UPZIM(I) END DO DO I = 1, NZERO UPZIM(I) = ARRLOK(IARR(I)) END DO DO I = 1, NZERO ARRLOK(I) = XYONEMIN(I) END DO DO I = 1, NZERO XYONEMIN(I) = ARRLOK(IARR(I)) END DO DO I = 1, NZERO ARRLOK(I) = XYONEMAX(I) END DO DO I = 1, NZERO XYONEMAX(I) = ARRLOK(IARR(I)) END DO DO I = 1, NZERO IARRLOK(I) = IZDOUBT(I) END DO DO I = 1, NZERO IZDOUBT(I) = IARRLOK(IARR(I)) END DO DO I = 1, NZERO IARRLOK(I) = ISHEET(I) END DO DO I = 1, NZERO ISHEET(I) = IARRLOK(IARR(I)) END DO WRITE(1,*) WRITE(1,*) NZERO,' ! nzero; now up; nzleft=',NZLEFT, $ ' ! "isheet" rel "icuttyp=2"' WRITE(1,*) DO I = 1, NZERO IF (IZDOUBT(I).EQ.0) THEN WRITE(1,*) UPZRE(I),UPZIM(I),ISHEET(I) ELSE IF (IZDOUBT(I).EQ.1) THEN WRITE(1,*) UPZRE(I),UPZIM(I),ISHEET(I),' %1' ELSE WRITE(1,*) UPZRE(I),UPZIM(I),ISHEET(I),' %2' END IF END DO WRITE(1,*) WRITE(1,*) NZERO,' ! nzero; now "xyonemin,xyonemax"' WRITE(1,*) DO I = 1, NZERO IF (IZDOUBT(I).EQ.0) THEN WRITE(1,*) XYONEMIN(I),XYONEMAX(I) ELSE IF (IZDOUBT(I).EQ.1) THEN WRITE(1,*) XYONEMIN(I),XYONEMAX(I),' %1' ELSE WRITE(1,*) XYONEMIN(I),XYONEMAX(I),' %2' END IF END DO WRITE(1,*) WRITE(1,*) NZERO,' ! nzero; now "dlndu,derivdu"' WRITE(1,*) DO I = 1, NZERO ISHEETTYP = ISHEET(I) UP = DCMPLX(UPZRE(I),UPZIM(I)) CALL GRDERIVDU ( UP, SLASK,CSLASK, DLNDU,DERIVDU ) IF (IZDOUBT(I).EQ.0) THEN WRITE(1,*) DLNDU,DERIVDU ELSE IF (IZDOUBT(I).EQ.1) THEN WRITE(1,*) DLNDU,DERIVDU,' %1' ELSE WRITE(1,*) DLNDU,DERIVDU,' %2' END IF END DO WRITE(1,*) WRITE(1,*) NZERO,' ! nzero; now "group-u"' WRITE(1,*) DO I = 1, NZERO ISHEETTYP = ISHEET(I) UP = DCMPLX(UPZRE(I),UPZIM(I)) CALL GRDERIV ( UP, SLASK,CSLASK, DLNDO,DERIVDO, DLNDU,DERIVDU ) CSLASK = U0*UP + OMEG*(- DEXP(DLNDO-DLNDU)*DERIVDO/DERIVDU) IF (IZDOUBT(I).EQ.0) THEN WRITE(1,*) CSLASK ELSE IF (IZDOUBT(I).EQ.1) THEN WRITE(1,*) CSLASK,' %1' ELSE WRITE(1,*) CSLASK,' %2' END IF END DO IF (NRUS0.GT.0) THEN C extrapolation for the zeros WRITE(11,*) WRITE(11,*) ' --------------------------------------- ' WRITE(11,*) ' --------------------------------------- ' WRITE(11,*) ' --------------------------------------- ' WRITE(11,*) DO I = 1, NZERO IF (ISHEET(I).EQ.0) GOTO 101 ! actually impossible here IXROW = 1 ZC(IAZ_TAB(1,1)) = DCMPLX(UPZRE(I),UPZIM(I)) DO WHILE (.TRUE.) IXROW = IXROW + 1 IF (IIFSMSOME.NE.0) THEN IF (EPSFSUB.LE.0.0D0) HSTEPFSUB = HSTEPFSUB0 / MHQ(IXROW) END IF IF (IXROW.GT.1) THEN MHQEFF = MHQNOW * MHQ(IXROW) CALL RMODEL ( 0,0, MHQEFF ) END IF ISHEETTYP = ISHEET(I) XYSLASK = DMAX1 ( XYONEMIN(I)+0.2D0*XYONENORM , XYONENORM ) XYSLASK = DMIN1(XYONEMAX(I),XYSLASK) NZEROQ = 0 ICORNTEST = 0 UPIQQQ1 = UPZIM(I) - XYSLASK UPIQQQ2 = UPZIM(I) + XYSLASK CALL ZEROFIND (CDETSUB,CDETSUBCIRC, $ UPZRE(I)-XYSLASK,UPZRE(I)+XYSLASK,UPIQQQ1,UPIQQQ2, $ DDSUPMIN,DDSUPMAX,DDSUPFAKTMAX, DERIVDRELMAX, $ DRELMAX,DRELGOALFAKT, ARGDIFSTORE, SLIMLEFT, $ DLOGREDUCFAKT,DDELT,XYTOL, IOPTION,ICHECK,11, $ nzerodim, NZEROQ,NZLEFTQ, UPZREQ,UPZIMQ, $ XYONEMINQ,XYONEMAXQ, NELMAX,NPLINMAX,IASFCALL, $ IDZHFCALL,ICHKFCALL,INELCOMPR,IASBACK,IDZHBACK, $ ICHKBACK) IF ( NZEROQ.NE.1 .OR. NZLEFTQ.NE.0 ) THEN ISHEET(I) = 0 GOTO 101 END IF ZC(1) = DCMPLX ( UPZREQ(1), UPZIMQ(1) ) C zc(iaz_tabnew(1)) = zc(1) ! done since iadr00_tabnew = 0 JCOLX = MIN0(IXROW,JRUS0) IF (IXROW.GE.IRUS0) THEN OKSUB = .TRUE. ELSE OKSUB = .FALSE. END IF CALL XPOLBUL ( 1,1, IXTYP0, 1.0D0,EPS0, $ HQI(1,IXROW),FIHQI,QHQI, 2,JCOLX, $ IAZ_TABNEW,IAZ_VALTERM,IAZ_TAB,IAZ_YTAB, $ ZC,OKSUB ) IF ( OKSUB .OR. IXROW.GE.NRUS0 ) THEN ZC(1) = ZC(IAZ_VALTERM(1)) UPZRE(I) = DREAL(ZC(1)) UPZIM(I) = DIMAG(ZC(1)) GOTO 101 END IF END DO 101 CONTINUE END DO WRITE(1,*) WRITE(1,*) WRITE(1,*) NZERO,' ! nzero; extrapolated zeros' WRITE(1,*) DO I = 1, NZERO IF (ISHEET(I).NE.0) THEN WRITE(1,*) UPZRE(I),UPZIM(I) ELSE WRITE(1,*) ' % failure' END IF END DO DO I = 1, 6, 5 ! write on file & terminal WRITE(I,*) WRITE(I,*) WRITE(I,*)' INCLUDING EXTRAPOLATION' WRITE(I,*) ' nelmax,nplinmax = ',NELMAX,NPLINMAX WRITE(I,*) ' iasfcall,idzhfcall,ichkfcall = ',IASFCALL, $ IDZHFCALL,ICHKFCALL WRITE(I,*) ' inelcompr = ',INELCOMPR WRITE(I,*) ' iasback,idzhback,ichkback = ',IASBACK, $ IDZHBACK,ICHKBACK END DO END IF DO I = 1, 6, 5 ! write on file & terminal WRITE(I,*) WRITE(I,*) WRITE(I,*) ' nbig was ',NBIG WRITE(I,*) ' nzero,nzleft = ',NZERO,NZLEFT WRITE(I,*) ' nelmax,nplinmax = ',NELMAX,NPLINMAX WRITE(I,*) ' iasfcall,idzhfcall,ichkfcall = ',IASFCALL, $ IDZHFCALL,ICHKFCALL WRITE(I,*) ' inelcompr = ',INELCOMPR WRITE(I,*) ' iasback,idzhback,ichkback = ',IASBACK,IDZHBACK, $ ICHKBACK IF (RHO(NOL).GT.0.0D0) THEN WRITE(I,*) WRITE(I,*) ' upbpa was ',UPBPA IF (ISOLID(NOL).NE.0) THEN WRITE(I,*) ' upbpb was ',UPBPB END IF IF (ICUTTYP0.EQ.1) THEN WRITE(I,*) ' nzdoubt1,nzdoubt2 ',NZDOUBT1,NZDOUBT2 WRITE(I,*) ' (pertaining to "xyonemin regions" ', $ 'around axis-hyperbolic branch-cuts)' END IF END IF IF (RHO(0).GT.0.0D0) THEN WRITE(I,*) WRITE(I,*) ' upbpasurf was ',UPBPASURF IF (ISOLID(0).NE.0) THEN WRITE(I,*) ' upbpbsurf was ',UPBPBSURF END IF IF (ICUTTYP0.NE.2) STOP 341552 ! not implemented END IF WRITE(I,*) END DO CLOSE(1) CLOSE(11) IF (IOPTION.EQ.2) CLOSE(99) END DO STOP END INTEGER FUNCTION ISHEETLOK (IPS,ISHEETTYP0) C for isheettyp IMPLICIT NONE INTEGER IPS,ISHEETTYP0 IF (IPS.EQ.1) THEN IF (ISHEETTYP0.EQ.1) THEN ISHEETLOK = 2 ELSE IF (ISHEETTYP0.EQ.2) THEN ISHEETLOK = 1 ELSE IF (ISHEETTYP0.EQ.3) THEN ISHEETLOK = 4 ELSE IF (ISHEETTYP0.EQ.4) THEN ISHEETLOK = 3 END IF ELSE IF (IPS.EQ.2) THEN IF (ISHEETTYP0.EQ.1) THEN ISHEETLOK = 3 ELSE IF (ISHEETTYP0.EQ.3) THEN ISHEETLOK = 1 ELSE IF (ISHEETTYP0.EQ.2) THEN ISHEETLOK = 4 ELSE IF (ISHEETTYP0.EQ.4) THEN ISHEETLOK = 2 END IF ELSE IF (IPS.EQ.3) THEN IF (ISHEETTYP0.EQ.1) THEN ISHEETLOK = 4 ELSE IF (ISHEETTYP0.EQ.4) THEN ISHEETLOK = 1 ELSE IF (ISHEETTYP0.EQ.2) THEN ISHEETLOK = 3 ELSE IF (ISHEETTYP0.EQ.3) THEN ISHEETLOK = 2 END IF ELSE ISHEETLOK = ISHEETTYP0 END IF RETURN END INTEGER FUNCTION IAZ_TABNEW (IHR) IMPLICIT NONE INTEGER IHR INCLUDE 'IADR00Z.INC' IAZ_TABNEW = IADR00_TABNEW + IHR RETURN END INTEGER FUNCTION IAZ_VALTERM (IHR) IMPLICIT NONE INTEGER IHR INCLUDE 'IADR00Z.INC' IAZ_VALTERM = IADR00_VALTERM + IHR RETURN END INTEGER FUNCTION IAZ_TAB (IHR,J) IMPLICIT NONE INTEGER IHR,J INCLUDE 'IADR00Z.INC' IAZ_TAB = IADR00_TAB + (J-1)*1 + IHR RETURN END INTEGER FUNCTION IAZ_YTAB (IHR,J) IMPLICIT NONE INTEGER IHR,J INCLUDE 'IADR00Z.INC' IAZ_YTAB = IADR00_YTAB + (J-1)*1 + IHR RETURN END SUBROUTINE CORNTESTSUB (X1,X2,Y1,Y2, SKIPRECT) IMPLICIT NONE COMMON /CORNTEST/ ICORNTEST, XYHYPA,XYHYPB INTEGER ICORNTEST REAL*8 XYHYPA,XYHYPB REAL*8 X1,X2,Y1,Y2 LOGICAL SKIPRECT REAL*8 SLASK, SLASKMIN,SLASKMAX IF (ICORNTEST.EQ.0) THEN SKIPRECT = .FALSE. RETURN ELSE GOTO (1,2,3,4,5,6,7) ICORNTEST STOP 4516 1 IF (X1.GE.0.0D0) THEN IF ( X1*Y1 .GT. XYHYPA ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF ELSE IF (X2.GT.0.0D0) STOP 341 IF ( X2*Y2 .GT. XYHYPA ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF END IF RETURN 2 IF (X1.GE.0.0D0) THEN IF ( X2*Y2 .LT. XYHYPA ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF ELSE IF (X2.GT.0.0D0) STOP 341 IF ( X1*Y1 .LT. XYHYPA ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF END IF RETURN 3 IF (X1.GE.0.0D0) THEN IF ( X1*Y1 .GT. XYHYPB ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF ELSE IF (X2.GT.0.0D0) STOP 341 IF ( X2*Y2 .GT. XYHYPB ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF END IF RETURN 4 IF (X1.GE.0.0D0) THEN IF ( X2*Y2 .LT. XYHYPB ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF ELSE IF (X2.GT.0.0D0) STOP 341 IF ( X1*Y1 .LT. XYHYPB ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF END IF RETURN 5 SLASK = DMIN1(XYHYPA,XYHYPB) IF (X1.GE.0.0D0) THEN IF ( X1*Y1 .GT. SLASK ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF ELSE IF (X2.GT.0.0D0) STOP 341 IF ( X2*Y2 .GT. SLASK ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF END IF RETURN 6 SLASK = DMAX1(XYHYPA,XYHYPB) IF (X1.GE.0.0D0) THEN IF ( X2*Y2 .LT. SLASK ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF ELSE IF (X2.GT.0.0D0) STOP 341 IF ( X1*Y1 .LT. SLASK ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF END IF RETURN 7 SLASKMIN = DMIN1(XYHYPA,XYHYPB) SLASKMAX = DMAX1(XYHYPA,XYHYPB) IF (X1.GE.0.0D0) THEN SLASK = X2*Y2 IF ( SLASK .LT. SLASKMIN ) THEN SKIPRECT = .TRUE. ELSE SLASK = X1*Y1 IF ( SLASK .GT. SLASKMAX ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF END IF ELSE IF (X2.GT.0.0D0) STOP 341 SLASK = X1*Y1 IF ( SLASK .LT. SLASKMIN ) THEN SKIPRECT = .TRUE. ELSE SLASK = X2*Y2 IF ( SLASK .GT. SLASKMAX ) THEN SKIPRECT = .TRUE. ELSE SKIPRECT = .FALSE. END IF END IF END IF RETURN END IF END