C*****************************************************************

      PROGRAM RINIT

C*****************************************************************
      
      INCLUDE 'rinit_inc.f'
      
      TIME=0.D0
      TIME0=TIME
      DSIG(0)=0.

      FH3=1.

C FH3 = A CONSTANT THAT CAN BE USED TO MODIFY SHEAR

C*****************************************************************

C     VARSINAISET INPUT ARVOT

      
      READ(*,*) DUMMYLINE
      WRITE(*,*) ' LKM,ETA,SEED'
      READ(*,*) LKM,ETA,ISEED
      WRITE(*,*) ' KAPPA/OMEGA, WZ/OMEGA'
      READ(*,*) KAPPAFAC,WZFAC
      WRITE(*,*) 'ROO'
      READ(*,*) ROO

      W0=DSQRT(GM/ETA/ETA/ETA)
      OMEGA=W0
      W2=2.D0*OMEGA
      KAPPA=KAPPAFAC*OMEGA
      OMEGA2=OMEGA*OMEGA
      KAPPA2=KAPPA*KAPPA
      WZ=W0*WZFAC
      WZ2=WZ*WZ
      PERIOD=2.D0*PII/OMEGA

      READ(*,*) xl,yl
      WRITE(*,*) xl,yl
      READ(*,*) cink,ceks
      WRITE(*,*) cink,ceks


C     INITIALIZE NAG-RANDOM NUMBER GENERATOR
      DUMMY=RAN1(ISEED)
      IIDUM=0


C*********************************************************************
C     MASS DISTRIBUTION
C*********************************************************************

      WRITE(*,*) 'DISCRETE (1) OR CONTINUOUS DISTRIBUTION (2)'
      READ(*,*) VAL

C     ****************************************
C     DISCRETE
C     ****************************************

      IF(VAL.EQ.1) THEN
         MJAK=1
         WRITE(*,*) 'NUMBER OF SIZE GROUPS'
         READ(*,*) NSIG
         WRITE(*,*) ' SIG, NUMBER, CIFAC, CEFAC'
         IND=0
         DO 5 I=1,NSIG
            READ(*,*) SIG,NUM(I),CIFAC(I),CEFAC(I)
            IF(I.EQ.1) SIGA=SIG
            DSIG(I)=1.0001*SIG
            DO 6 J=1,NUM(I)
               IND=IND+1
               S(IND)=SIG
               XMASS(IND)=4.*PII/3.*ROO*S(IND)**3
               IS(IND)=I
               MYY(IND)=XMASS(IND)
               IF(ROO.LE.0) 
     +              MYY(IND)=4.*PII/3.*900.*S(IND)**3
 6          CONTINUE
 5       CONTINUE
         WRITE(*,*) ' NTOT ',IND
         LKM=IND
      ENDIF

C     ****************************************
C     CONTINUOUS DISTRIBUTION
C     ****************************************

      IF(VAL.EQ.2) THEN
         MJAK=2
         WRITE(*,*) ' NUMBER OF SIZE-BINS'
         READ(*,*) NSIG
         WRITE(*,*) ' Q,SIGA,SIGB'
         READ(*,*) Q,SIGA,SIGB
         QQ=Q
         VAK= DLOG(SIGB/SIGA)/NSIG
         DO 10 I=1,NSIG
            DSIG(I)=DEXP(VAK*I)*SIGA*1.00001
 10      CONTINUE
         EXPO=1-Q
         S1=SIGA**EXPO
         S2=SIGB**EXPO
         DO 19 I=1,LKM
            SS=(1.*(I-1)/(LKM-1)*(S2-S1)+S1)**(1./EXPO)
            MYY(I)=(SS/SIGA)**3
            S(I)=SS
            XMASS(I)=4.*PII/3.*ROO*S(I)**3
            MYY(I)=XMASS(I)
            IF(ROO.LE.0) 
     +           MYY(I)=4.*PII/3.*900.*S(I)**3

 19      CONTINUE
         WRITE(*,*) 'SMIN,SMAX',S(1),S(LKM)
      ENDIF

      WRITE(*,*) 'DSIG'
      WRITE(*,*) (DSIG(I),I=0,NSIG)


c---------------------------------------------
c Optical depths
c---------------------------------------------


C     ******* IF XL=0 USE YL AS TAU


      cxmax=xl

      IF(XL.EQ.0) THEN
         PA=0.
         DO 4567 I=1,LKM
            PA=PA+PII*S(I)*S(I)
 4567    CONTINUE
         XL=SQRT(PA/YL/4.)
         YL=XL
         CYMAX=0.9999*YL
         CX1(1)=-CYMAX
         CX2(1)=CYMAX
         cxmax=cymax
         
         WRITE(6,*) 'XL,YL -> '
         WRITE(6,*) 'TAU=',PA/4./XL/YL
         WRITE(6,*) 'XL=YL=',XL
         WRITE(6,*) 'CXMAX=CYMAX=',CYMAX
      ENDIF


      TAU0=0.
      NTEST=0
      MTOT_SIZE(0)=0.
      DO 29 I=1,NSIG
         NUM(I)=0
         TAU(I)=0.
         MKES(I)=0.
         MTOT_SIZE(I)=0.
 29   CONTINUE
      DO 30 I=1,LKM
         IS(I)=NSIG
         DO 31 J=1,NSIG
            IF((S(I).LE.DSIG(J)).AND.(S(I).GT.DSIG(J-1)))
     +           THEN
               IS(I)=J
               NUM(J)=NUM(J)+1
               TAU(J)=TAU(J)+PII*S(I)*S(I)/CXMAX/YL/4.
               TAU0=TAU0+PII*S(I)*S(I)/CXMAX/YL/4.
               MKES(J)=MKES(J)+MYY(I)
               MTOT_SIZE(J)=MTOT_SIZE(J)+MYY(I)
               MTOT_SIZE(0)=MTOT_SIZE(0)+MYY(I)
               GO TO 30
            ENDIF
 31      CONTINUE
 30   CONTINUE




      FACTEMP=KAPPA2-4.D0*OMEGA2
      SHEAR=-FACTEMP/2.D0/OMEGA
      AVAK=2.D0*XL
      CVAK=2.D0*YL
      VVAK=FACTEMP/OMEGA*XL
      BVAK=VVAK*TIME


      WRITE(*,*) ' simple_180112-SIMULATION     '
      WRITE(*,*) ' DISTANCE FROM PLANET (10E3 KM)',ETA/1.D6
      WRITE(*,*) ' KAPPA, WZ',KAPPAFAC,WZFAC
      WRITE(*,*) ' ORIGINAL     CXMAX CYMAX',CXMAX,CYMAX
      WRITE(*,*) ' ORIGINAL     INK2 EKS2',CINK,CEKS
      WRITE(*,*) ' BOX-SIZE     XL YL    ',XL,YL
      WRITE(*,*) '  '
      WRITE(*,*) 'GROUP   SIG1-SIG2   MA-MB   NN   MKES   TAU'
      WRITE(*,*) '  '

      DO 40 I=1,NSIG
         SIG1=DSIG(I-1)
         SIG2=DSIG(I)
         IF(NUM(I).NE.0) MKES(I)=MKES(I)/NUM(I)
         MA=(SIG1/SIGA)**3
         MB=(SIG2/SIGA)**3
         WRITE(*,41) I,SIG1,SIG2,MA,MB,NUM(I),MKES(I),TAU(I)
 40   CONTINUE
 41   FORMAT(I6,2F8.3,2F10.3,I8,F10.3,F12.6)





C**************************************************************
C     CREATE RADIUS AND VELOCITY VECTORS
C**************************************************************

      I=0
      lkmr=1

      LASK=0
      INK2=0.D0
      EKS2=0.D0
      CXMI=-xl
      CXMA=xl

      DO 111 IR=1,lkm/2
         I=I+2
         CIFAC0=1.
         CEFAC0=1.
         IF(VAL.EQ.1) THEN
            CIFAC0=CIFAC(IS(I))
            CEFAC0=CEFAC(IS(I))
         ENDIF
         
 112     CONTINUE

C     ******* CREATE PARTICLES IN PAIRS SO THAT
C     UVAK=SUMM (MYY XP) = 0
C     WVAK=SUMM (MYY(YP+1.5OMEGAX)) = 0
C     THIS IS TRUE IF FOR EACH PAIR
C     MYY1 EKS1 = MYY2 EKS2
C     DELTA1 = DELTA2 + PII
C     RAYLEIGH DISTRIBUTION OF EKS,INC (SEE TRULSEN 1972)

C     FIRST PARTICLE OF THE PAIR INDEX=I
         SRAN=RAN1(IIDUM)
         CX=CXMI+SRAN*(CXMA-CXMI)
         SRAN=RAN1(IIDUM)
         CY=-CYMAX+2.*SRAN*CYMAX
         SRAN=RAN1(IIDUM)
         DELTA=2.D0*PII*SRAN
         SRAN=RAN1(IIDUM)
         GAMMA=2.D0*PII*SRAN
         SRAN=RAN1(IIDUM)
         INK=CIFAC0*CINK*DSQRT(DLOG(1.D0/(1.D0-SRAN)))
         SRAN=RAN1(IIDUM)
         EKS=CEFAC0*CEKS*DSQRT(DLOG(1.D0/(1.D0-SRAN)))
         R(I,1)=CX+EKS*DCOS(DELTA)
         R(I,2)=CY-2.D0*OMEGA/KAPPA*EKS*DSIN(DELTA)
         R(I,3)=INK*DSIN(GAMMA)
         RP(I,1)=-KAPPA*EKS*DSIN(DELTA)
         RP(I,2)=-2.D0*W0*EKS*DCOS(DELTA)+
     +        FH3*CX*(KAPPA2-4.D0*OMEGA2)/2.D0/OMEGA
         RP(I,3)=INK*WZ*DCOS(GAMMA)

C     SECOND PARTICLE OF THE PAIR INDEX=I-1

         SRAN=RAN1(IIDUM)
         CXB=CXMI+SRAN*(CXMA-CXMI)
         SRAN=RAN1(IIDUM)
         CYB=-CYMAX+2.*SRAN*CYMAX
         EKSB=MYY(I)/MYY(I-1)*EKS
         INKB=MYY(I)/MYY(I-1)*INK
         
         R(I-1,1)=CXB+EKSB*DCOS(DELTA+PII)
         R(I-1,2)=CYB-2.D0*OMEGA/KAPPA*EKSB*DSIN(DELTA+PII)
         R(I-1,3)=INKB*DSIN(GAMMA+PII)
         RP(I-1,1)=-KAPPA*EKSB*DSIN(DELTA+PII)
         RP(I-1,2)=-2.D0*W0*EKSB*DCOS(DELTA+PII)+
     +        FH3*CXB*(KAPPA2-4.D0*OMEGA2)/2.D0/OMEGA
         RP(I-1,3)=INKB*WZ*DCOS(GAMMA+PII)
         

C     **************************************
C     MOVE PAIR INTO CALCULATION BOX

         N=0
         M=0
         IF(R(I,1).GE.0) N=-INT(R(I,1)/AVAK+0.5)
         IF(R(I,1).LT.0) N=INT(-R(I,1)/AVAK+0.5)
         R(I,1)=R(I,1)+N*AVAK
         IF(R(I,2).GE.0) M=-INT(R(I,2)/CVAK+0.5)
         IF(R(I,2).LT.0) M=INT(-R(I,2)/CVAK+0.5)
         R(I,2)=R(I,2)+M*CVAK
         RP(I,2)=RP(I,2)+N*VVAK
         
         N=0
         M=0
         IF(R(I-1,1).GE.0) N=-INT(R(I-1,1)/AVAK+0.5)
         IF(R(I-1,1).LT.0) N=INT(-R(I-1,1)/AVAK+0.5)
         R(I-1,1)=R(I-1,1)+N*AVAK
         IF(R(I-1,2).GE.0) M=-INT(R(I-1,2)/CVAK+0.5)
         IF(R(I-1,2).LT.0) M=INT(-R(I-1,2)/CVAK+0.5)
         R(I-1,2)=R(I-1,2)+M*CVAK
         RP(I-1,2)=RP(I-1,2)+N*VVAK

         LASK=LASK+1
         IF(LASK.GT.lkm) THEN
            LASK=1
            CINK=CINK*1.2
            WRITE(*,*) 'CINK',CINK
         ENDIF

C     **************************************
C     CHECK FOR OVERLAPPING PARTICLES
C     (ALSO OVER BORDERS)

         DO 800 J=1,I-1
            EROX=R(I,1)-R(J,1)
            EROY=R(I,2)-R(J,2)
            EROZ=R(I,3)-R(J,3)
            IF(DABS(EROY).GT.YL) EROY=2.D0*YL-DABS(EROY)
            IF(DABS(EROX).GT.XL) EROX=2.D0*XL-DABS(EROX)
            ERO=EROX*EROX+EROY*EROY+EROZ*EROZ
            SKER=(S(I)+S(J))**2
            IF(ERO.LT.SKER) GOTO 112
            IF(J.EQ.I-1) GOTO 800
            EROX=R(I-1,1)-R(J,1)
            EROY=R(I-1,2)-R(J,2)
            EROZ=R(I-1,3)-R(J,3)
            IF(DABS(EROY).GT.YL) EROY=2.D0*YL-DABS(EROY)
            IF(DABS(EROX).GT.XL) EROX=2.D0*XL-DABS(EROX)
            ERO=EROX*EROX+EROY*EROY+EROZ*EROZ
            SKER=(S(I-1)+S(J))**2
            IF(ERO.LT.SKER) GOTO 112
 800     CONTINUE
         LASK=LASK-1

C     *******  PAIR WAS ACCEPTED

         INK2=INK2+INK*INK+INKB*INKB
         EKS2=EKS2+EKS*EKS+EKSB*EKSB
         XPRO(I)=CX
         YPRO(I)=CY
         XPRO(I-1)=CXB
         YPRO(I-1)=CYB
         IF(MOD(I,50).EQ.0) WRITE(*,*) I
 111  CONTINUE
      
C     ******* OUTPUT MEAN VALUES FOR EACH RADIAL GROUP
      
      WRITE(*,*) 'SQRT(INK2) =',SQRT(INK2/lkm)
      WRITE(*,*) 'SQRT(EKS2) =',SQRT(EKS2/lkm)


C*****************************************************************
C     POSITIONS AND VELOCITIES CALCULATED
C     CHECK THAT THE SYSTEM IS AT REST
C*****************************************************************


C*****************************************************************
C     CREATE HIGHER DENSITIES FOR 2D-CASE
C     IF CINK < 0 -> PUT VERTICAL POSITIONS TO ZERO
C     IMPACTS TAKE CARE OF OVERLAPPS (HOPEFULLY)

      IDIME=3
      IF(CINK.LE.0) IDIME=2
      IF(CINK.LT.0) THEN
         DO 2222 I=1,LKM
            R(I,3)=0.
            RP(I,3)=0.
 2222    CONTINUE
      ENDIF
      
      
C***********************************************
C     17/01/03 ASSIGN SPINS

      DO 1234 I=1,LKM
         ROTA(I,1)=0.
         ROTA(I,2)=0.
         ROTA(I,3)=QZINIT*S(I)*W0
 1234 CONTINUE

C***********************************************

      UVAK=0.D0
      WVAK=0.D0
      ZVAK=0.D0
      ZPVAK=0.D0
      MSUM=0.
      EPOT=0.
      EKIN=0.

      DO 102 I=1,LKM
         UVAK=UVAK+RP(I,1)*MYY(I)
         WVAK=WVAK+(RP(I,2)+SHEAR*R(I,1))*MYY(I)
         ZVAK=ZVAK+R(I,3)*MYY(I)
         ZPVAK=ZPVAK+RP(I,3)*MYY(I)
         MSUM=MSUM+MYY(I)
         EKIN=EKIN+0.5*MYY(I)*(RP(I,1)**2+RP(I,2)**2+RP(I,3)**2)
         EPOT=EPOT+0.5*MYY(I)*(FACTEMP*R(I,1)**2+WZ2*R(I,3)**2)
 102  CONTINUE
      UVAK=UVAK/W0/XL/MSUM
      WVAK=UVAK/W0/XL/MSUM
      ZVAK=ZVAK/MSUM
      ZPVAK=ZPVAK/W0/XL/MSUM
      UVAK_ALL=UVAK
      WVAK_ALL=WVAK
      ZVAK_ALL=ZVAK
      ZPVAK_ALL=ZPVAK
      ETOT=EKIN+EPOT
      ETOT0=ETOT

      WRITE(*,*) 'UVAK JA WVAK'
      WRITE(*,*) UVAK,WVAK
      WRITE(*,*) 'ZVAK JA ZPVAK'
      WRITE(*,*) ZVAK,ZPVAK
      WRITE(*,*) 'EKIN JA EPOT'
      WRITE(*,*) EKIN,EPOT

      
      TMASS=0.D0
      DO 13 I=1,LKM
         TMASS=TMASS+XMASS(I)
 13   CONTINUE
      WRITE(*,*) 'RINIT: TMASS',TMASS


C---------------------------------------------------
C     WRITE FINPOSNEW-FILE (FOR SIMPLE)
C---------------------------------------------------
      
      OPEN(UNIT=50,FILE='iniposnew',STATUS='UNKNOWN')
      WRITE(6,*) 'OUTPUT POSITIONS TO FINPOSNEW'
      WRITE(6,*) 'TIME,LKM',TIME,LKM
      WRITE(50,*) TIME ,LKM
      WRITE(50,*) (R(I,1),R(I,2),R(I,3),I=1,LKM)
      WRITE(50,*) (RP(I,1),RP(I,2),RP(I,3),I=1,LKM)
      WRITE(50,*) (MYY(I),I=1,LKM)
      WRITE(50,*) (S(I),I=1,LKM)
      WRITE(50,*) (XMASS(I),I=1,LKM)
      WRITE(50,*) (ROTA(I,1),ROTA(I,2),ROTA(I,3),I=1,LKM)
      WRITE(50,*) (IS(I),I=1,LKM)
      WRITE(50,*) ETA
      WRITE(50,*) KAPPAFAC,WZFAC
      WRITE(50,*) XL,YL
      CLOSE(50)

C---------------------------------------------------

      END PROGRAM



C*****************************************************

      FUNCTION RAN1(IDUM)

C*****************************************************
      DOUBLE PRECISION RAN1
      DIMENSION R(97)
      PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6)
      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6)
      PARAMETER (M3=243000,IA3=4561,IC3=51349)
C     DATA IFF /0/
      SAVE
C     IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
      IF (IDUM.LT.0) THEN
       WRITE(*,*) 'IDUM-SEED',IDUM
        IFF=1
        IX1=MOD(IC1-IDUM,M1)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IX1,M2)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX3=MOD(IX1,M3)
        DO 11 J=1,97
          IX1=MOD(IA1*IX1+IC1,M1)
          IX2=MOD(IA2*IX2+IC2,M2)
          R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
11      CONTINUE
        IDUM=1
      ENDIF
      IX1=MOD(IA1*IX1+IC1,M1)
      IX2=MOD(IA2*IX2+IC2,M2)
      IX3=MOD(IA3*IX3+IC3,M3)
      J=1+(97*IX3)/M3
C      IF(J.GT.97.OR.J.LT.1) PAUSE
      RAN1=R(J)
      R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1

C      WRITE(99,*) 'RND1: ',IDUM,RAN1
      RETURN
      END

