SUBROUTINE M700FM
* (XL,YL,SIDEX,SIDEY,NITER,MR,MI,XCMPND,CYCLE,
* RATE,ROOT,COL,MCPAX,MCPAY)
C
C COPYRIGHT (C) 1985-1992 BY ART MATRIX. ALL RIGHTS RESERVED.
C PROGRAMMER BC LAST REVISION 2/16/90
C * A PETER PROGRAM - F
C
C***********************************************************************
C***********************************************************************
C
INCLUDE 'FSIGNAL.INC'
INCLUDE 'RC.INC'
REAL*8 XL,YL,SIDEX,SIDEY,AR,AI,MR,MI,XCMPND
REAL*8 EPS
REAL*8 X,X1,X2,X3,XLAST,HOLDX
REAL*8 Y,Y1,Y2,Y3,YLAST
REAL*8 FCMP,FXA,FXB,FXC,FYA,FYB,FYC
REAL*8 GCMP,GXA,GXB,GXC,GYA,GYB,GYC
REAL*8 F,FX,FY,G,GX,GY
REAL*8 DETJAC,CHANGE,SPHERE,BALL
C
INTEGER*2 RATE(0:MCPAX,0:MCPAY)
INTEGER*2 ROOT(0:MCPAX,0:MCPAY)
INTEGER*2 COL (0:MCPAX,0:MCPAY)
INTEGER*4 ROOTX,ROOTY,CYCLE
INTEGER*4 IX,IY,NITER,CMPND,COUNT1,COUNT2,MOSTIT
C
LOGICAL*4 HEADRM
C
PARAMETER (MAXINT = 32767)
C
C***********************************************************************
C STATEMENT FUNCTIONS :
C***********************************************************************
C
F() = -2*X3 + 5*X1*Y2 - X1 + AR
FX() = -6*X2 + 5*Y2 - 1
FY() = 10*X1*Y1
C
G() = 2*Y3 - 5*X2*Y1 - Y1 + AI
GX() = -10*X1*Y1
GY() = 6*Y2 - 5*X2 - 1
C
C***********************************************************************
C***********************************************************************
C
C TRAP AND IGNORE OVERFLOWS.
CALL XOFLOW
C
C***********************************************************************
C CONSTANT ASSINGNMENTS :
C***********************************************************************
C
C -MINIMUM RADIUS BETWEEN POINTS CONSIDERED CONVERGENT.
SPHERE = 1.D-4
C
C -MAKE SURE SPHERE IS LESS THAN 1/2 A PIXEL IN RADIUS.
BALL = (MIN(SIDEX,SIDEY) / MCPAX) / 2
IF (SPHERE .GT. BALL) THEN
SPHERE = BALL
END IF
C
C EVILQ POSTER!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C SPHERE=.000053129D0
C COMPUTED FROM SIDEX = .38656
C MCPAX = 3638
C CROWN OF THORNS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SPHERE=8.5779105D-7
C COMPUTED FROM SIDEX = .20586
C MCPAX = 1200
C
C***********************************************************************
C***********************************************************************
C
C
C - TOO SMALL TO DIVIDE BY :
EPS = 1.D-5
C
C -NUMBER OF TIMES FUNCTION IS COMPOUNDED.
C CMPND = (INPUT VIA PARAMETER LINE ARGUEMENT.)
C
CMPND = XCMPND
WRITE (6,*) 'CMPND ',CMPND
C
C -INITIALIZE ROOT FREQUENCY ARRAY.
DO IY = 0,MCPAY,1
DO IX = 0,MCPAX,1
ROOT(IX,IY) = 0
END DO
END DO
C
C***********************************************************************
C MAIN OUTER LOOPS. LOOP THROUGH EACH POINT IN GRID (IX,IY)
C***********************************************************************
C
CALL NOERR
C
C DO 2000 IY = 177,MCPAY,1
DO 2000 IY = 0,MCPAY,1
C
AI = YL + (IY*SIDEY)/MCPAY
C
CALL SHOWIY (10,IY,AI)
IF (RC .EQ. 1) RETURN
C
C DO 1000 IX = 161,161,1
DO 1000 IX = 0,MCPAX,1
C
C WRITE (6,*) 'IY,IX',IY,IX
C
AR = XL + (IX*SIDEX)/MCPAX
C
C***********************************************************************
C NEWTONS METHOD CALCULATIONS FOR EACH POINT IN IX,IY GRID.
C***********************************************************************
C -ASSIGNMENTS REFRESHED PRIOR TO EACH POINT'S CALCULATIONS
C
XLAST = 1E10
YLAST = 1E10
X = MR
Y = MI
C
C -ITERATION LOOP; # OF TIMES TO APPLY NEWTON'S METHOD BEFORE GIVING UP.
C
DO 400 COUNT1 = 1,NITER
C
X1 = X
Y1 = Y
X2 = X1*X1
X3 = X2*X1
Y2 = Y1*Y1
Y3 = Y2*Y1
C
FCMP= F()
GCMP= G()
FXB = FX()
FYB = FY()
GXB = GX()
GYB = GY()
C
DO 200 COUNT2 = 2,CMPND
C
X1 = FCMP
Y1 = GCMP
X2 = X1*X1
X3 = X2*X1
Y2 = Y1*Y1
Y3 = Y2*Y1
C
FCMP= F()
GCMP= G()
FXA = FX()
FYA = FY()
GXA = GX()
GYA = GY()
C
FXC = FXB * FXA + FYB * GXA
GXC = GXB * FXA + GYB * GXA
FYC = FXB * FYA + FYB * GYA
GYC = GXB * FYA + GYB * GYA
C
FXB = FXC
FYB = FYC
GXB = GXC
GYB = GYC
C
200 CONTINUE
C
FCMP = FCMP - X
GCMP = GCMP - Y
C
C - M711,712,713,711A,712A DO NOT USE THE 4 FOLLOWING CORRECT LINES !
C - 2/16/90
FXB = FXB - 1
C M721,2,3 FYB = FYB - 1
C M721,2,3 GXB = GXB - 1
GYB = GYB - 1
C
DETJAC = (FXB*GYB - GXB*FYB)
C WRITE (6,*) FXB*GYB,GXB*FYB
C WRITE (6,*) 'DETJAC',DETJAC
IF (ABS(DETJAC) .LT. EPS) GOTO 600
C
HOLDX = X - (GYB*FCMP - FYB*GCMP) / DETJAC
Y = Y - (FXB*GCMP - GXB*FCMP) / DETJAC
X = HOLDX
CHANGE = ABS(XLAST - X) + ABS(YLAST - Y)
IF (CHANGE .LE. SPHERE) GOTO 800
XLAST = X
YLAST = Y
400 CONTINUE
C
RATE(IX,IY) = 0
GOTO 1000
C
600 RATE(IX,IY) = -1
C WRITE (6,*) 'GOTTA -1'
GOTO 1000
C
800 RATE(IX,IY) = COUNT1
ROOTX = (X - XL)*MCPAX/SIDEX
ROOTY = (Y - YL)*MCPAY/SIDEY
C
IF (ROOTX .LT. 0 .OR. ROOTX .GT. MCPAX .OR.
* ROOTY .LT. 0 .OR. ROOTY .GT. MCPAY) GOTO 1000
C
ROOT(ROOTX,ROOTY) = ROOT(ROOTX,ROOTY) + 1
C
GOTO 1000
C
C***********************************************************************
C***********************************************************************
C
1000 CONTINUE
2000 CONTINUE
C
C***********************************************************************
C COMBINE ROOT AND RATE DATA.
C***********************************************************************
C
MOSTIT = -40000
DO IY = 0,MCPAY,1
DO IX = 0,MCPAX,1
IF (RATE(IX,IY) .GT. MOSTIT) THEN
MOSTIT = RATE(IX,IY)
IF (MOSTIT .GE. NITER) GOTO 3000
END IF
END DO
END DO
C
3000 CONTINUE
C
HEADRM = .TRUE.
CYCLE = INT(MAXINT/MAX(1,MOSTIT))
C
DO IY = 0,MCPAY,1
DO IX = 0,MCPAX,1
IF (ROOT(IX,IY) .GE. CYCLE) THEN
COL(IX,IY) = CYCLE * RATE(IX,IY) + (CYCLE - 1)
HEADRM = .FALSE.
ELSE
COL(IX,IY) = CYCLE * RATE(IX,IY) + ROOT(IX,IY)
END IF
END DO
END DO
C
WRITE (6,*) 'CYCLE SIZE & MOST ITERATIONS : ',CYCLE,MOSTIT
C
IF (HEADRM) THEN
WRITE (6,*) 'ENOUGH HEADROOM EXISTED.'
ELSE
WRITE (6,*) 'NOT ENOUGH HEADROOM.'
END IF
C
C***********************************************************************
C***********************************************************************
C
CALL DRVID ('FM ',102,1,2,CYCLE)
C
RETURN
END