C @(#)chaos.for	17.1.1.1 (ESO-DMD) 01/25/02 17:14:56
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
      PROGRAM CHAOS
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  program CHAOS			version 1.00	901031
C  K. Banse				ESO - Garching
C
C.KEYWORDS
C  dynamic systems, bulk data frame
C
C.PURPOSE
C  create a MIDAS image
C
C.ALGORITHM
C  use MIDAS interfaces
C
C.INPUT/OUTPUT
C  the following keywords are used:
C
C IN_A/C/1/60                  	name of new image file
C DEFAULT/C/1/1         	= Y, all image info given in command line
C                  		= N, use frame given in key P2 for that
C P3/C/1/60                  	info frame, if DEFAULT = N
C INPUTI/I/1/7                 	NAXIS + NPIX
C INPUTR/R/1/12                	START + STEP
C P4/C/1/20                  	type of data function
C                  		= LORENZ, ELFLY, JULIA, SYMM
C P5/C/1/60                  	coefficients for data function
C
C.VERSIONS
C  1.00		from version 1.20 of crfram.for
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      NAXIS,CUNLEN
      INTEGER*8    PNTR
      INTEGER      IAV,IDUM,N,NCOEF,NO,STAT
      INTEGER      INO,IBUF(7),NPIX(6)
      INTEGER      NLO,RNO,UNIT(1),MADRID(1)
      INTEGER      NPIX1,NPIX2
C
      REAL         CUTS(4),COEFS(10),FMIN,FMAX,RDUM
C
      DOUBLE PRECISION START(6),STEP(6),DBUF(12)
      DOUBLE PRECISION DCOEFS(3),DDUM
C
      CHARACTER    CUNIT*112,IDENT*72,NAME*60,INFRAM*60
      CHARACTER    CBUF*60,DEFAUL*1
      CHARACTER    FUNCY*4
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      COMMON   /VMR/ MADRID
C
      DATA      CUNIT /'none given '/, CBUF /' '/
      DATA      NPIX      /6*1/,       START /6*0.D0/, STEP /6*1.D0/
      DATA      NAME      /' '/,       INFRAM /' '/
C
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C  get into MIDAS
      CALL STSPRO('CHAOS')
C
C  get name and type of new frame and DEFAULT flag
      CALL STKRDC('IN_A',1,1,60,IAV,NAME,UNIT,NLO,STAT)
      CALL STKRDC('DEFAULT',1,1,1,IAV,DEFAUL,UNIT,NLO,STAT)
      CALL UPCAS(DEFAUL,DEFAUL)
C
C  image info from command line or reference frame
      IF (DEFAUL.NE.'Y') THEN
         CALL STKRDC('P3',1,1,60,IAV,INFRAM,UNIT,NLO,STAT)  !get ref. frame
C
C  read all relevant descriptors
         CALL STFOPN(INFRAM,D_OLD_FORMAT,0,F_IMA_TYPE,RNO,STAT)
         CALL STDRDI(RNO,'NAXIS',1,1,IAV,NAXIS,UNIT,NLO,STAT)
         CALL STDRDI(RNO,'NPIX',1,NAXIS,IAV,NPIX,UNIT,NLO,STAT)
         CALL STDRDD(RNO,'START',1,NAXIS,IAV,START,UNIT,NLO,STAT)
         CALL STDRDD(RNO,'STEP',1,NAXIS,IAV,STEP,UNIT,NLO,STAT)
         CALL STDRDC(RNO,'IDENT',1,1,72,IAV,IDENT,UNIT,NLO,STAT)
         CALL STDRDC(RNO,'CUNIT',1,1,112,IAV,CUNIT,UNIT,NLO,STAT)
      ELSE
C
C     get NAXIS, NPIX, START + STEP (max 6 dimensions)
         CALL STKRDI('INPUTI',1,7,IAV,IBUF,UNIT,NLO,STAT)
         CALL STKRDD('INPUTD',1,12,IAV,DBUF,UNIT,NLO,STAT)
C
C     set descriptor values
         NAXIS = IBUF(1)
         IF ((NAXIS.LT.1).OR.(NAXIS.GT.6))
     +      CALL STETER(1,'invalid NAXIS ...')      !invalid NAXIS ...
C
         DO 800 N=1,NAXIS
            NPIX(N) = IBUF(N+1)
            START(N) = DBUF(N)
            STEP(N) = DBUF(N+NAXIS)
800      CONTINUE
      ENDIF
C
C  finally look for coefficients + function type
      CALL STKRDC('P4',1,1,4,IAV,FUNCY,UNIT,NLO,STAT)
      CALL UPCAS(FUNCY,FUNCY)
      CALL STKRDC('P5',1,1,60,IAV,CBUF,UNIT,NLO,STAT)
C
      IDENT = 'artificial image '
      CUNLEN = (NAXIS+1)*16
C
C  create new frame
      CALL STIPUT(NAME,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXIS,NPIX,START,STEP,IDENT,
     +            CUNIT(1:CUNLEN),PNTR,INO,STAT)
      NO = NPIX(1) * NPIX(2)
      CALL NL(MADRID(PNTR),NO)                   !set everything to 0.0
C
C  branch according to function type
C 
C
      NPIX1  =  NPIX(1)
      NPIX2  =  NPIX(2)
C  Lorenz attractors
      IF (FUNCY(1:2).EQ.'LO') THEN
         CALL GENCNV(CBUF,2,7,IDUM,COEFS,DDUM,NCOEF) !x,y,z + h,a,b + no_iter
         IF (NCOEF.LT.6) 
     +      CALL STETER(1,'we need 6 coefficients...')
         IF (NCOEF.LT.7) COEFS(7) = 4000.
         CALL LOR(MADRID(PNTR),NPIX1,NPIX2,COEFS,FMIN,FMAX,
     +                                              START,STEP)   
C 
C  chaotic dynamic system  - Electronic Fly)
      ELSE IF (FUNCY(1:3).EQ.'ELF') THEN
         CALL GENCNV(CBUF,2,6,IDUM,COEFS,DDUM,NCOEF)    !a,b,c,d,e + no_iter
         IF (NCOEF.LT.5) 
     +      CALL STETER(1,'we need 5 coefficients...')
         IF (NCOEF.LT.6) COEFS(6) = 1000000.
         CALL ELFLY(MADRID(PNTR),NPIX1,NPIX2,COEFS,FMIN,FMAX,
     +                                               START,STEP)  
C 
C  Mandelbrot set
      ELSE IF (FUNCY(1:1).EQ.'M') THEN
         CALL GENCNV(CBUF,4,3,IDUM,RDUM,DCOEFS,NCOEF)    !x,y + no_iter
         IF (NCOEF.LT.2) 
     +      CALL STETER(1,'we need 2 coefficients...')
         IF (NCOEF.LT.3) DCOEFS(3) = 100.
         CALL MANDEL(MADRID(PNTR),NPIX1,NPIX2,START,STEP,DCOEFS,
     +               FMIN,FMAX)
C 
C  Julia set
      ELSE IF (FUNCY(1:1).EQ.'J') THEN
         CALL GENCNV(CBUF,4,3,IDUM,RDUM,DCOEFS,NCOEF)    !x,y + no_iter
         IF (NCOEF.LT.2) 
     +      CALL STETER(1,'we need 2 coefficients...')
         IF (NCOEF.LT.3) DCOEFS(3) = 100.
         CALL JULIA(MADRID(PNTR),NPIX1,NPIX2,START,STEP,DCOEFS,
     +              FMIN,FMAX)
C 
C  Colour wheel
      ELSE IF (FUNCY(1:3).EQ.'COL') THEN
         CALL COLWH(MADRID(PNTR),NPIX1,NPIX2,FMIN,FMAX)
C 
C  Symmetric chaos
      ELSE IF (FUNCY(1:3).EQ.'SYM') THEN
         CALL GENCNV(CBUF,2,7,IDUM,COEFS,DDUM,NCOEF)            !a,b,c,d
         CALL SYM(MADRID(PNTR),NPIX1,NPIX2,COEFS,FMIN,FMAX)
C 
      ELSE
         CALL STETER(9,'invalid function...')
      ENDIF
C
C  update cuts
      CUTS(1) = 0.
      CUTS(2) = 0.
      CUTS(3) = FMIN
      CUTS(4) = FMAX
      CALL STDWRR(INO,'LHCUTS',CUTS,1,4,UNIT,STAT)
C     
      IF (FUNCY(1:2).EQ.'LO') THEN
         CUNIT(1:) = 'Lorenz: '//CBUF
      ELSE IF (FUNCY(1:3).EQ.'ELF') THEN
         CUNIT(1:) = 'ElFly: '//CBUF
      ELSE IF (FUNCY(1:1).EQ.'M') THEN
         CUNIT(1:) = 'Mandelbrot set: '//CBUF
      ELSE IF (FUNCY(1:1).EQ.'J') THEN
         CUNIT(1:) = 'Julia set: '//CBUF
      ELSE IF (FUNCY(1:3).EQ.'COL') THEN
         CUNIT(1:) = 'Colour wheel '
      ELSE IF (FUNCY(1:3).EQ.'SYM') THEN
         CUNIT(1:) = 'Symmetric chaos: '//CBUF
      ENDIF
      CALL STDWRC(INO,'HISTORY',1,CUNIT,1,80,UNIT,STAT)
C
C  that's it folks...
      CALL STSEPI
      END

      SUBROUTINE LOR(A,NPIX1,NPIX2,RR,FMI,FMA,START,STEP)
C 
      IMPLICIT NONE
C 
      INTEGER   NPIX1,NPIX2
      INTEGER   ITLIM,N,NX,NY
      DOUBLE PRECISION  START(2), STEP(2)
C 
      REAL      A(NPIX1,NPIX2),RR(*),FMI,FMA
      REAL      X,Y,Z,XNEW,YNEW,ZNEW,H,AAAA,FRAC
      REAL      XXMIN,XXMAX,YYMIN,YYMAX
      REAL      XINC,YINC,XXX,YYY
C 
      XXMIN = START(1)
      XXMAX = START(1)+NPIX1*STEP(1)
      YYMIN = START(2)
      YYMAX = START(2)+NPIX2*STEP(2)
      XINC  = 1./STEP(1)
      YINC  = 1./STEP(2)
C 
      X = RR(1)
      Y = RR(2)
      Z = RR(3)
      H = RR(4)
      AAAA = RR(5)
      FRAC = RR(6)
      ITLIM = INT(RR(7)+0.5)
C 
      DO 1000 N=1,ITLIM
         XNEW = X + H*10.*(Y-X)
         YNEW = Y + H*((-X*Z)+AAAA*X-Y)
         ZNEW = Z + H*(X*Y-FRAC*Z)
         X = XNEW
         Y = YNEW
         Z = ZNEW
         IF ( (X.GE.XXMIN) .AND. (X.LE.XXMAX) .AND.
     +        (Y.GE.YYMIN) .AND. (Y.LE.YYMAX) ) THEN
            XXX = (X - XXMIN) * XINC
            YYY = (Y - YYMIN) * YINC
            NX = INT(XXX+0.5)
            NY = INT(YYY+0.5)
            A(NX,NY) = 1.0
         ENDIF
1000  CONTINUE
C 
      FMI = 0.
      FMA = 1.
C 
      RETURN
      END

      SUBROUTINE NL(A,NSIZE)
C 
      IMPLICIT NONE
C 
      INTEGER  NSIZE,N
C 
      REAL  A(NSIZE)
C 
      DO 500 N=1,NSIZE
         A(N) = 0.
500   CONTINUE
C 
      RETURN
      END

      SUBROUTINE ELFLY(A,NPIX1,NPIX2,RR,FMI,FMA,START,STEP)
C 
      IMPLICIT NONE
C 
      INTEGER   NPIX2,NPIX1
      INTEGER   ITLIM,N,NX,NY
      DOUBLE PRECISION START(2), STEP(2)
C 
      REAL      A(NPIX1,NPIX2),RR(*),FMI,FMA
      REAL      AA,BB,CC,DD,EE
      REAL      XXMIN,XXMAX,YYMIN,YYMAX
      REAL      XINC,YINC,X,Y,Z,X1,Y1,Z1,XXX,YYY
C 
      AA = RR(1)
      BB = RR(2)
      CC = RR(3)
      DD = RR(4)
      EE = RR(5)
      ITLIM = INT(RR(6)+0.5)
C 
      XXMIN = START(1)
      XXMAX = START(1)+NPIX1*STEP(1)
      YYMIN = START(2)
      YYMAX = START(2)+NPIX2*STEP(2)
      XINC  = 1./STEP(1)
      YINC  = 1./STEP(2)
      X = 0.                            !start wit (0.,0.,0.)
      Y = 0.
      Z = 0.
C 
      DO 1000 N=1,ITLIM
         X1 = SIN(AA*Y) - Z*COS(BB*X)
         Y1 = Z*SIN(CC*X) - COS(DD*Y)
         Z1 = EE*SIN(X)
         X = X1
         Y = Y1
         Z = Z1
         IF ( (X.GE.XXMIN) .AND. (X.LE.XXMAX) .AND.
     +        (Y.GE.YYMIN) .AND. (Y.LE.YYMAX) ) THEN
            XXX = (X - XXMIN) * XINC
            YYY = (Y - YYMIN) * YINC
            NX = INT(XXX+0.5)
            NY = INT(YYY+0.5)
            A(NX,NY) = A(NX,NY) + 1.0
C            IF (A(NX,NY).GT.FMA) FMA = A(NX,NY)
         ENDIF
1000  CONTINUE
C 
      FMI = 0.
      FMA = 10.
C 
      RETURN
      END

      SUBROUTINE MANDEL(OUT,NPIX1,NPIX2,START,STEP,DCOF,FMI,FMA)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.ALGORITHM
C  Mandelbrot set in the tranformation    z <- z*z + c
C  with z and c complex numbers
C  input is used as initial value for z, current (x,y) coord -> c
C 
C-------------------------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER NPIX1,NPIX2,IX,IY
      INTEGER K,KUP
C
      REAL OUT(NPIX1,NPIX2)  
      REAL FMI,FMA
C
      DOUBLE PRECISION START(2)  
      DOUBLE PRECISION STEP(2)  
      DOUBLE PRECISION DCOF(3)
      DOUBLE PRECISION UPPER  
C  
      DOUBLE PRECISION DK,X1,Y1,X2,Y2,X,Y,DX,DY,R
C
      FMI = 0.
      FMA = 0.
      KUP = 256
      DX = STEP(1)
      DY = STEP(2)
      UPPER = DCOF(3)
C 
C  loop over all lines
      Y = START(2)
      DO 400 IY=1,NPIX2
         X = START(1)
C 
C  loop over all pixels of a line
         DO 300 IX=1,NPIX1
            X1 = DCOF(1)                      !Re(c) -> Re(z)
            Y1 = DCOF(2)                      !Im(c) -> Im(z)
C 
C  iterate over max 256 pixel levels
            DO 100 K=1,KUP
               X2 = (X1*X1) - (Y1*Y1) + X          !Re(z~) = Re(z*z) + Re(c)
               Y2 = 2.*(X1*Y1) + Y                 !Im(z~) = Im(z*z) + Im(c)
               R = X2*X2 + Y2*Y2                   !Abs(z~)
               IF (R.GT.UPPER) THEN
                  DK = K - 1
                  OUT(IX,IY) = DK
                  IF (DK.GT.FMA) FMA = DK
                  GO TO 200
               ENDIF
               X1 = X2
               Y1 = Y2
100         CONTINUE                !if we come out here, leave pixel at 0.0
C 
200         X = X + DX
300      CONTINUE
C 
         Y = Y + DY
400      CONTINUE
C 
      RETURN
      END


      SUBROUTINE SYM(OUT,NPIX1,NPIX2,COF,FMI,FMA)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.ALGORITHM
C
C-------------------------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER NPIX1,NPIX2,IX,IY
      INTEGER N,LIMIT
C
      REAL OUT(NPIX1,NPIX2)
      REAL FMI,FMA
      REAL COF(7)
C 
      DOUBLE PRECISION A,B,C,D,RFY,RFX,Z
      DOUBLE PRECISION X,Y,XX,YY,XXX,YYY,XN,YN
C
      X = COF(1)
      Y = COF(2)
      A = COF(3)
      B = COF(4)
      C = COF(5)
      D = COF(6)
      LIMIT = COF(7)
      RFX = NPIX1 - 1.
      RFY = NPIX2 - 1.
C
      DO 500, N=1,LIMIT
         XX = X * X
         XXX = XX * X
         YY = Y * Y
         YYY = YY * Y
         Z = B * (XX + YY)
         XN = A*X + B*X*(XX + YY) + C*(XXX-3.0*X*YY) + D*(XX-YY)
         YN = A*Y + B*Y*(XX + YY) + C*(3.0*XX*Y-YYY) - 2.0*D*X*Y
         X = XN
         Y = YN
         IX = (X*RFX) + 1
         IF (IX.LT.1) THEN
            IX = 1
         ELSE IF (IX.GT.NPIX1) THEN
            IX = NPIX1
         ENDIF
         IY = (Y*RFY) + 1
         IF (IY.LT.1) THEN
            IY = 1
         ELSE IF (IY.GT.NPIX2) THEN
            IY = NPIX2
         ENDIF
         OUT(IX,IY) = OUT(IX,IY) + 1.
500   CONTINUE
C 
      FMI = OUT(1,1)
      FMA = FMI
      DO 800, IY=1,NPIX2
         DO 700, IX=1,NPIX1
            IF (OUT(IX,IY).LT.FMI) THEN
               FMI = OUT(IX,IY)
            ELSE IF (OUT(IX,IY).GT.FMA) THEN
               FMA = OUT(IX,IY)
            ENDIF
700      CONTINUE
800   CONTINUE
C 
      RETURN
      END


      SUBROUTINE JULIA(OUT,NPIX1,NPIX2,START,STEP,DCOF,FMI,FMA)
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.ALGORITHM
C  Julia set in the tranformation    z <- z*z + c
C  with z and c complex numbers
C  input -> c, current (x,y) coord are used as initial value for z
C 
C-------------------------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER NPIX1,NPIX2,IX,IY
      INTEGER K,KUP
C
      REAL OUT(NPIX1,NPIX2)  
      REAL FMI,FMA
C
      DOUBLE PRECISION START(2)  
      DOUBLE PRECISION STEP(2)  
      DOUBLE PRECISION DCOF(3)
      DOUBLE PRECISION UPPER  
C  
      DOUBLE PRECISION DK,X1,Y1,X2,Y2,X,Y,DX,DY,R
C
      FMI = 0.
      FMA = 0.
      KUP = 256
      DX = STEP(1)
      DY = STEP(2)
      UPPER = DCOF(3)
C 
C  loop over all lines
      Y = START(2)
      DO 400 IY=1,NPIX2
         X = START(1)
C 
C  loop over all pixels of a line
         DO 300 IX=1,NPIX1
            X1 = X                      !xcoord -> Re(z)
            Y1 = Y                      !ycoord -> Im(z)
C 
C  iterate over max 256 pixel levels
            DO 100 K=1,KUP
               X2 = (X1*X1) - (Y1*Y1) + DCOF(1)       !Re(z~) = Re(z*z) + Re(c)
               Y2 = 2.*(X1*Y1) + DCOF(2)              !Im(z~) = Im(z*z) + Im(c)
               R = X2*X2 + Y2*Y2                   !Abs(z~)
               IF (R.GT.UPPER) THEN
                  DK = K - 1
                  OUT(IX,IY) = DK
                  IF (DK.GT.FMA) FMA = DK
                  GO TO 200
               ENDIF
               X1 = X2
               Y1 = Y2
100         CONTINUE                !if we come out here, leave pixel at 0.0
C 
200         X = X + DX
300      CONTINUE
C 
         Y = Y + DY
400      CONTINUE
C 
      RETURN
      END

      SUBROUTINE COLWH(A,NPIX1,NPIX2,FMIN,FMAX)
C 
      IMPLICIT NONE
C 
      INTEGER   NPIX1,NPIX2
      INTEGER   NX,NY,KX,KY
      INTEGER   MAXRAD,KOFF
C 
      REAL      PHI,R,X,Y
      REAL      A(1),FMIN,FMAX
C 
      FMIN = -3.2
      FMAX = 3.142
      KX = NPIX1/2
      KY = NPIX2/2
      MAXRAD = NINT(KX * 0.80)
      MAXRAD = MAXRAD * MAXRAD               !take square of radius
      KOFF = 0
C 
      DO 5000 NY=1,NPIX2
         Y = NY - KY
C 
         DO 4000 NX=1,NPIX1
            X = NX - KX
            R = X*X + Y*Y
            IF (R .GT. MAXRAD) THEN
               PHI = FMIN
            ELSE
               PHI = ATAN2(X,Y)
            ENDIF
            A(KOFF+NX) = PHI
4000     CONTINUE
         KOFF = KOFF + NPIX1
5000  CONTINUE
C 
      RETURN
      END