COMMON /RNDM/IU,IV, /LIMIT/XL,YL, /RESULT/NO,NRES,FMAX
      DIMENSION X(200),Y(200),AV(2),XL(2),YL(2),RSLT(50,5)
      CHARACTER *8 FIN,FOUT
      WRITE(*,*)'PROGRAM TO FIT A LOGARITHMIC SPIRAL TO EMPIRICAL DATA'
      WRITE(*,*)'WHEN ORIGIN HAS BEEN DISPLACED BY Cx AND Cy'
      WRITE(*,*)'ALGORITHM AND PROGRAM BY SK MISHRA, DEPT. OF ECONOMICS'
      WRITE(*,*)'NORTH-EASTERN HILL UNIVERSITY, SHILLONG (INDIA)'
      WRITE(*,*) '*****************************************************'
      WRITE(*,*)' '
      WRITE(*,*)'NAME THE INPUT FILE AND OUTPUT FILE'
      READ(*,*)FIN,FOUT
C     NI AND NO ARE THE UNIT NOS. OF THE INPUT AND OUTPUT FILES
      NI=7
      NO=8
      FMAX=0
      NRES=0
      OPEN(NI, FILE=FIN)
      READ(NI,*)N
      DO 1 I=1,N
      READ(NI,*) X(I),Y(I)
    1 CONTINUE
      CLOSE(NI)
      WRITE(*,*)'FEED RANDOM NUMBER SEED (a 4-digit odd number) '
      READ(*,*) IU
      WRITE(*,*)'FROM GRAPH OF THE SPIRAL GUESS NARROWEST RANGE Cx AND'
      WRITE(*,*)'Cy IN WHICH CENTRE OF THE SPIRAL LIES - JUST GUESS'
      WRITE(*,*)'LOWER AND UPPER LIMITS OF CX ?'
      READ(*,*)XL(1),XL(2)
      WRITE(*,*)'LOWER AND UPPER LIMITS OF CY ?'
      READ(*,*)YL(1),YL(2)
      AVX=0.0
      AVY=0.0
      DO 2 IM=1,N
      AVX=AVX+X(IM)
      AVY=AVY+Y(IM)
    2 CONTINUE
      AV(1)=AVX/N
      AV(2)=AVY/N
C     CALL RANDOMIZED BOX ALGORITHM FOR OPTIMIZATION
      CALL BOX(N,X,Y,AV,FOUT,FOPT,CX,CY,XLO,XHI,YLO,YHI)
      WRITE(*,*)'-----------------------------------------------------'
      WRITE(*,*)'TERMINAL VALUE OF FUNCTION = ',FOPT
      WRITE(*,*)'TERMINAL CX AND CY VALUES'
      WRITE(*,*) 'CX=',CX,'   CY=',CY
      WRITE(*,*)'RANGE : CX(',XLO,XHI,');   CY(',YLO,YHI,')'
      WRITE(*,*)'RESULTS OF OPTIMALITY SEARCH STORED IN FILE=',FOUT
C     OPTIMALITY SEARCH ENDED.
      OPEN(NI,FILE=FOUT)
C      TO BE READ FROM FILE FOUT : A(1),A(2),F,ALPHA,BETA
      DO 3 I=1,NRES
      READ(NI,*)(RSLT(I,J),J=1,5)
    3 CONTINUE
      CLOSE(NI)
      DO 4 I=1,NRES-1
      DO 4 II=I,NRES
      IF(RSLT(I,1).GT.RSLT(II,1)) THEN
      DO 5 J=1,5
      TEMP=RSLT(I,J)
      RSLT(I,J)=RSLT(II,J)
      RSLT(II,J)=TEMP
    5 CONTINUE
      ENDIF
    4 CONTINUE
      WRITE(*,*)' '
      WRITE(*,*)'******************************************************'
      WRITE(*,*)'RESULTS OF SEARCH IN ORDER OF IMPROVING R-SQUARE'
      WRITE(*,*)'-----------------------------------------------------'
      WRITE(*,*)'      R-SQUARE           Cx             Cy
     *a              b'
      DO 6 I=1,NRES
      WRITE(*,7)(RSLT(I,J),J=1,5)
    6 CONTINUE
      WRITE(*,*)'      R-SQUARE           Cx             Cy
     *a              b'
      WRITE(*,*)'---------------------- END     ------------------ '
      WRITE(*,*)'PROGRAM TERMINATED SUCCESSFULLY'
    7 FORMAT(1X,5F15.8)
      END

      SUBROUTINE BOX(NN,XX,YY,AV,FOUT,FOPT,CX,CY,XLO,XHI,YLO,YHI)
C     ADAPTED FROM THE CODES DEVELOPED BY JA RICHARDSON, ARIZONA STATE
C     UNIV. TEMPLE, ARIZONA LISTED IN J. KUESTER & J. MIZE (1973)
C     OPTIMIZATION TECHNIQUES WITH FORTRAN, MCGRAW HILLS INC. PP.368-385
      COMMON /RNDM/IU,IV, /LIMIT/XL,YL /RESULT/NO,NRES,FMAX
      DIMENSION X(30,2),R(30,2),F(30),G(2),H(2),XC(2),AV(2),XL(2),YL(2)
      DIMENSION XX(200),YY(200)
      CHARACTER *8 FOUT
      INTEGER GAMMA
C     WRITE(*,*)'FEED N,M,K,ITMAX,IC,IPRINT'
C     N=2 = NO. OF VARIABLES;
C     K=> N+1 : RANDOM NUMBERS GENERATED WITHIN THE RANGE OF VARIABLES
C     M = NO OF CONSTRAINTS
      N=2
      M=2
      K=30
C     ITMAX IS NO. OF MAX ITERATIONS
      WRITE(*,*)'NO. OF ITERATIONS TO PERFORM (SAY 10 TO 30) '
      READ(*,*) ITMAX
      WRITE(*,*)'-------------- COMPUTING PLEASE WAIT --------'
C     IPRINT=0 MEANS NO INTERMEDIATE PRINT REQUIRED.
      IPRINT=1
      VALU=0
      ALPHA=1.3
      BETA=0.00001
      GAMMA=5
      DELTA=0.00001
C     WRITE(*,*)'FEED (X(1,J),J=1,N), THE INITIAL GUESSES=AVERAGES ?'
      CALL CONST(G,H)
      DO 4 J=1,N
      IF((AV(J).GE.H(J)).OR.(AV(J).LE.G(J))) THEN
      AV(J)=(H(J)+G(J))/2.0
      ENDIF
      X(1,J)=AV(J)
    4 CONTINUE
C     OPEN THE OUTPUT FILE
      OPEN(NO,FILE=FOUT)
C     ITERATIVE RANDOMIZED SEARCH BEGINS
C     ITERATIVE LOOP FOR SHRINKAGE OF RANGE OF XLO,XHI,YLO AND YHI
      DO 5 ITL=1,30
C     ITERATIVE LOOP FOR OPTIMIZATION
      DO 6 ITER=1,50
      DO 1 J=1,N
      DO 1 I=2,K
      CALL RANDOM(RAND)
      R(I,J)=RAND
    1 CONTINUE
      CALL CONSX(N,M,K,ITMAX,ALPHA,BETA,GAMMA,DELTA,
     *X,R,F,IT,IEV1,IEV2,G,H,XC,IPRINT,XX,YY,NN)
      VL=F(1)
      DO 2 JJ=2,K
      IF(F(JJ).GT.VL) THEN
      JMAX=JJ
      VL=F(JJ)
      ENDIF
    2 CONTINUE
      IF(VALU.LT.VL) THEN
      VALU=VL
      DO 3 J=1,N
      X(1,J)=X(JMAX,J)
    3 CONTINUE
      F(1)=VL
      ENDIF
    6 CONTINUE
C     LOCAL OPTIMUM IS FOUND
      FOPT=F(1)
      CX=X(1,1)
      CY=X(1,2)
      XLO=XL(1)
      XHI=XL(2)
      YLO=YL(1)
      YHI=YL(2)
      WRITE(*,*)'CURRENT BEST VALUE OF FUNCTION = ',F(1)
      WRITE(*,*)'CURRENT -- CX AND CY VALUES'
      WRITE(*,*) 'CX=',X(1,1),'   CY=',X(1,2)
      WRITE(*,*)'IN THE RANGE:CX(',XL(1),XL(2),'); CY(',YL(1),YL(2),')'
C     SHRINKAGE OF RANGE OF DECISION VARIABLES : XLO,XHI,YLO AND YHI
      CALL RANDOM(RAND)
      XL(1)=XL(1)-(XL(1)-X(1,1))*2*RAND/3.0
      CALL RANDOM(RAND)
      XL(2)=XL(2)-(XL(2)-X(1,1))*2*RAND/3.0
      CALL RANDOM(RAND)
      YL(1)=YL(1)-(YL(1)-X(1,2))*2*RAND/3.0
      CALL RANDOM(RAND)
      YL(2)=YL(2)-(YL(2)-X(1,2))*2*RAND/3.0
    5 CONTINUE
      CLOSE(NO)
      RETURN
      END
C     ---------------------------------------------------
      SUBROUTINE CONSX(N,M,K,ITMAX,ALPHA,BETA,GAMMA,DELTA,X,R,F,IT,IEV1,
     *IEV2,G,H,XC,IPRINT,XX,YY,NN)
      DIMENSION X(30,2),R(30,2),F(30),G(2),H(2),XC(2),XX(200),YY(200)
      INTEGER GAMMA
C      WRITE(*,*)'ENTERED IN CONSX'
C      DO 3 II=1,NN
C      WRITE(*,*) II,XX(II),YY(II),  'IN CONSX'
C   3 CONTINUE
      IBREAK=0
      IT =1
      KODE=0
      IF(M-N) 20,20,10
   10 KODE=1
   20 CONTINUE
      DO 40 II=2,K
      DO 30 J=1,N
   30 X(II,J)=0.0
   40 CONTINUE
      DO 65 II=2,K
      DO 50 J=1,N
      I=II
      CALL CONST(G,H)
      X(II,J)=G(J)+R(II,J)*(H(J)-G(J))
   50 CONTINUE
      K1=II
      CALL CHECK(N,M,X,G,H,I,KODE,XC,DELTA,K1)
      IF(II-2) 51,51,55
   51 IF(IPRINT) 52,65,52
   52 CONTINUE
C      WRITE(*,*)'CO-ORDINATES OF COMPLEX'
C      I0=1
C      WRITE(*,*) (I0,J,X(I0,J),J=1,N)
   55 IF(IPRINT) 56,65,56
   56 CONTINUE
C      WRITE(*,*)(II,J,X(II,J),J=1,N)
   65 CONTINUE
      K1=K
      DO 70 I=1,K
      CALL FUNC(N,X,F,I,XX,YY,NN)

   70 CONTINUE
      KOUNT=1
      IF(IPRINT) 72,80,72
   72 CONTINUE
C      WRITE(*,*)(J,F(J),J=1,K)
   80 IEV1=1
      DO 100 ICM=2,K
      IF(F(IEV1)-F(ICM)) 100,100,90
   90 IEV1=ICM
  100 CONTINUE
C
      IEV2=1
      DO 120 ICM=2,K
      IF(F(IEV2)-F(ICM)) 110,110,120
  110 IEV2=ICM
  120 CONTINUE
C
      IF(F(IEV2)-(F(IEV1)+BETA)) 140,130,130
  130 KOUNT=1
      GOTO 150
  140 KOUNT=KOUNT+1
      IF(KOUNT-GAMMA) 150,240,240
  150 CALL CENTR(N,IEV1,XC,X,K1)
      DO 160 JJ=1,N
  160 X(IEV1,JJ) = (1.0+ALPHA)*XC(JJ)-ALPHA*X(IEV1,JJ)

      I=IEV1
      CALL CHECK(N,M,X,G,H,I,KODE,XC,DELTA,K1)
      CALL FUNC(N,X,F,I,XX,YY,NN)
C
  170 IEV2=1
      DO 190 ICM=2,K
      IF(F(IEV2)-F(ICM)) 190,190,180
  180 IEV2=ICM
  190 CONTINUE
      IF(IEV2-IEV1) 220,200,220
  200 DO 210 JJ=1,N
      X(IEV1,JJ)=(X(IEV1,JJ)+XC(JJ))/2.0
  210 CONTINUE
      I=IEV1
      CALL CHECK(N,M,X,G,H,I,KODE,XC,DELTA,K1)
      CALL FUNC(N,X,F,I,XX,YY,NN)
      IBREAK=IBREAK+1
      IF(IBREAK.GT.10) GOTO 220
      GOTO 170
  220 CONTINUE
      IF(IPRINT) 230,228,230
  230 CONTINUE

  228 IT=IT+1

      IPRINT=0
      IF(IT/5.0.EQ.INT(IT/5)) IPRINT=1
      IF(IT-ITMAX) 80,80,240
  240 RETURN
      END
C       ---------------------------------------------------
      SUBROUTINE FUNC(N,X,F,I,XX,YY,NN)
      DIMENSION X(30,2),F(30),XX(200),YY(200),A(2)
      DO 1 J=1,N
      A(J)=X(I,J)
    1 CONTINUE
      CALL FUNCT(NN,XX,YY,A,FF)
      F(I)=FF
      RETURN
      END
C     ---------------------------------------------------
      SUBROUTINE CENTR(N,IEV1,XC,X,K1)
C     FINDS THE CENTROID

      DIMENSION X(30,2),XC(2)
C      WRITE(*,*)'ENTERED IN CENTR'
      DO 20 J=1,N
      XC(J)=0.0
      DO 10 IL=1,K1
   10 XC(J)=XC(J)+X(IL,J)
      RK=K1
   20 XC(J)=(XC(J)-X(IEV1,J))/(RK-1.0)
      RETURN
      END
C     ------------------------------------------------
      SUBROUTINE CONST(G,H)
      COMMON /LIMIT/XL,YL
      DIMENSION G(2),H(2),XL(2),YL(2)
C     SPECIFY THE CONSTRAINTS
C      WRITE(*,*)'ENTERED IN CONST'
C      LOWER AND UPPER LIMITS OF PARAMETERS
      G(1)=XL(1)
      H(1)=XL(2)
      G(2)=YL(1)
      H(2)=YL(2)
      RETURN
      END
C     ----------------------------------------------------
      SUBROUTINE CHECK(N,M,X,G,H,I,KODE,XC,DELTA,K1)
      DIMENSION X(30,2),G(2),H(2),XC(2)
C      WRITE(*,*)'ENTERED IN CHECK'
   10 KT=0
      CALL CONST(G,H)
C     CHECK AGAINST EXPLICIT CONSTRAINTS
      DO 50 J=1,N
      IF(X(I,J)-G(J)) 20,20,30
   20 X(I,J)=G(J)+DELTA
      GOTO 50
   30 IF(H(J)-X(I,J)) 40,40,50
   40 X(I,J)=H(J)-DELTA
   50 CONTINUE
      IF(KODE) 110,110,60
C     CHECK AGAINST THE IMPLICIT FUNCTION
   60 NN=N+1
      IF(NN.LE.M) THEN
      DO 100 J=NN,M
      CALL CONST(G,H)
      IF(X(I,J)-G(J)) 80,70,70
   70 IF(H(J)-X(I,J)) 80,100,100
   80 IEV1=I
      KT=1
      CALL CENTR(N,IEV1,XC,X,K1)
      DO 90 JJ=1,N
      X(I,JJ)=(X(I,JJ)+XC(JJ))/2.0
   90 CONTINUE
  100 CONTINUE
      ENDIF
      IF(KT) 110,110,10
  110 RETURN
      END

      SUBROUTINE RANDOM(RAND)
      COMMON /RNDM/IU,IV
      IV=IU*65539
      IF(IV.LT.0) THEN
      IV=IV+2147483647+1
      ENDIF
      RAND=IV
      IU=IV
      RAND=RAND*0.4656613E-09
      RETURN
      END

      SUBROUTINE FUNCT(N,X,Y,A,F)
      COMMON /RESULT/NO,NRES,FMAX
      DIMENSION A(2),X(200),Y(200),R(200),Q(200),AK(200),ANG(200)
      DIMENSION D(2),C(2,2),XO(200),YO(200)
C      WRITE(*,*)'ENTERED IN FUNCT**'
C     EVALUATE THE FUNCTION

      PI=4.0*ATAN(1.0)
C     PRESERVING X AND Y
      DO 12 I=1,N
      XO(I)=X(I)
      YO(I)=Y(I)
   12 CONTINUE
      DO 1 I=1,N
      X(I)=X(I)-A(1)
      Y(I)=Y(I)-A(2)
      R(I)=SQRT(X(I)**2+Y(I)**2)
C      WRITE(*,*) X(I),Y(I),R(I)
    1 CONTINUE
      DO 6 I=1,N
      Q(I)=1
      IF((X(I).LT.0.0).AND.(Y(I).GE.0.0)) Q(I)=2
      IF((X(I).LT.0.0).AND.(Y(I).LT.0.0)) Q(I)=3
      IF((X(I).GE.0.0).AND.(Y(I).LT.0.0)) Q(I)=4
    6 CONTINUE
      DO 7 I=1,N-1
      DO 8 II=I+1,N
      IF(R(I).GT.R(II)) THEN
      TEMP=Q(I)
      Q(I)=Q(II)
      Q(II)=TEMP
      TEMP=X(I)
      X(I)=X(II)
      X(II)=TEMP
      TEMP=Y(I)
      Y(I)=Y(II)
      Y(II)=TEMP
      TEMP=R(I)
      R(I)=R(II)
      R(II)=TEMP
      ENDIF
    8 CONTINUE
    7 CONTINUE
      KK=0
      AK(1)=KK
      DO 5 I=2,N
      AK(I)=KK
      IF(Q(I).LT.Q(I-1)) THEN
      KK=KK+1
      AK(I)=KK
      ENDIF
    5 CONTINUE
      DO 9 I=1,N
      IF(ABS(X(I)).LE.1.0E-30) THEN
      IF(Y(I).GT.0) ANGO=PI/2.0
      IF(Y(I).LT.0) ANGO=3.0*PI/2.0
      ENDIF
      IF(ABS(X(I)).GT.1.0E-30) ANGO=ATAN(Y(I)/X(I))
      ANG(I)=INT(Q(I)/2)*PI+ANGO
C      ANG(I)=ATAN(Y(I)/X(I))
C      IF(Q(I).EQ.2) ANG(I)=PI+ANG(I)
C      IF(Q(I).EQ.3) ANG(I)=PI+ANG(I)
C      IF(Q(I).EQ.4) ANG(I)=2*PI+ANG(I)
      ANG(I)=ANG(I)+2.0*PI*AK(I)
    9 CONTINUE
C     LEAST SQUARES ESTIMATION
      SSR=0
      SR=0
      SANG=0
      SSANG=0
      SRANG=0
      DO 10 I=1,N
      RL=LOG(R(I))
      SR=SR+RL
      SANG=SANG+ANG(I)
      SSANG=SSANG+ANG(I)**2
      SRANG=SRANG+RL*ANG(I)
      SSR=SSR+RL**2
   10 CONTINUE
      D(1)=SR
      D(2)=SRANG
      C(1,1)=N
      C(1,2)=SANG
      C(2,1)=SANG
      C(2,2)=SSANG
C     MATRIX INVERSION OF C
      TEMP=C(1,1)
      C(1,1)=C(2,2)
      C(2,2)=TEMP
      C(1,2)=-C(1,2)
      C(2,1)=-C(2,1)
      DET=C(1,1)*C(2,2)-C(1,2)*C(2,1)
      C(1,1)=C(1,1)/DET
      C(1,2)=C(1,2)/DET
      C(2,1)=C(2,1)/DET
      C(2,2)=C(2,2)/DET
C     INVERTED MATRIX C-INVERSE IS POST-MULTIPLIED BY D
      AA=C(1,1)*D(1)+C(1,2)*D(2)
      BB=C(2,1)*D(1)+C(2,2)*D(2)
C     ROOT MEAN SQUARES OF ERRORS
      SMS=0.0
      DO 11 I=1,N
      SMS=SMS+(LOG(R(I))-AA-BB*ANG(I))**2
  11  CONTINUE
C      RMS=SQRT(SMS/(N))
      VARR=SSR/N-(SR/N)**2
      RSQUARE=1.0-(SMS/N)/VARR
C      TA=AA/(RMS*SQRT(C(1,1)))
C      TB=BB/(RMS*SQRT(C(2,2)))
C     TAKING ANTILOG(AA)
      AA= EXP(AA)
      F=RSQUARE
      ALPHA=AA
      BETA=BB
C      WRITE(*,*)'A,B,F',AA,BB,F
C     RESUMING ORIGINAL X
      IF(F.GT.FMAX) THEN
      FMAX=F
      NRES=NRES+1
      WRITE(NO,*)F,A(1),A(2),ALPHA,BETA
      ENDIF

      DO 13 I=1,N
      X(I)=XO(I)
      Y(I)=YO(I)
   13 CONTINUE
      RETURN
      END