*** This work is unpublished, please contact sealin2008@hotmail.com for more ******
*** detailed information. *********************************************************
***********************************************************************************
      PROGRAM MAIN
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(NX=512,PI=3.1415926535897932384626)
      DIMENSION X(0:NX),D(0:NX),W(0:NX),WX(0:2*NX-1),WY(0:2*NX-1)
      DIMENSION CT(0:2*NX-1),RCT(0:2*NX-1),TCT(0:2*NX-1),DEN(0:2*NX-1)

!     Initial D values
      RX = 0.35; TX = 1.0
      XL = 12.00; XIN = 0.85; XOUT = 4.0; DX = XL/NX
      XM = XL/2.; DT = 0.00001; NT = 4000
      DO I = 0, NX
C        X(I) = I*DX; D(I) = 0.5*SIN(2.*PI*X(I)/XL)
C        X(I) = I*DX; D(I) = 0.2*(X(I)-XM)*EXP(-(X(I)-XM)**2/2.)
        X(I) = I*DX; D(I) = 0.0
      ENDDO
      D(NX-1) = 0.1; D(NX-3) = - 0.1; D(NX-5) = 0.05; D(NX-7) = - 0.05 

      T = 0.0
C      OPEN(10,FILE='Dens.txt')
C      READ(10,*) T
C      DO I = 0, NX
C        READ(10,*) X(I),D(I)
C      ENDDO
C      READ(10,*) DT,NT
C      CLOSE(10)

      CALL WVAL(W,WX,WY,XIN,XOUT,DX,NX)
      CALL FFTVAL(CT,RCT,TCT,DEN,DX,NX,DT,RX,TX)
      ENG = ENERGY(D,W,WX,WY,RX,TX,DX,NX)

      OPEN(20,FILE='Evolution.txt')
      OPEN(30,FILE='Average.txt')

      WRITE(20,"(1X,A4,1X,E13.6,1X,A7,1x,E13.6)") 'Tao=',T,'Energy=',ENG
      SUM = 0.0
      DO J = 0, NX
        WRITE(20,'(2(1X,E13.6))') X(J),D(J)
        IF((J.EQ.0).OR.(J.EQ.NX)) THEN
          SUM = SUM + 0.5*ABS(D(J))
        ELSE
          SUM = SUM + ABS(D(J))
        ENDIF
      ENDDO
      SUM = SUM/NX
      WRITE(30,'(3(1X,E13.6))') T,SUM,ENG
      WRITE(*,*) T, SUM, ENG

      DO I = 1, 150
        T = T + DT*NT
        CALL STEP_FFT(D,WX,WY,DX,NX,NT,RX,TX,CT,RCT,TCT,DEN)
C        CALL STEP_FDM(D,W,WX,WY,DX,NX,DT,NT,RX,TX)
C        CALL STEP_DIR(D,W,WX,WY,DX,NX,DT,NT,RX,TX)
        ENG = ENERGY(D,W,WX,WY,RX,TX,DX,NX)
      WRITE(20,"(1X,A4,1X,E13.6,1X,A7,1x,E13.6)") 'Tao=',T,'Energy=',ENG
        SUM = 0.0
        DO J = 0, NX
          WRITE(20,'(2(1X,E13.6))') X(J),D(J)
          IF((J.EQ.0).OR.(J.EQ.NX)) THEN
            SUM = SUM + 0.5*ABS(D(J))
          ELSE
            SUM = SUM + ABS(D(J))
          ENDIF
        ENDDO
        SUM = SUM/NX
        WRITE(30,'(4(1X,E13.6))') T,SUM,ENG
        WRITE(*,'(4(1X,E13.6))') T, SUM, ENG
      ENDDO

      OPEN(10,FILE='Dens.txt')
      WRITE(10,'(1X,E13.6)') T
      DO I = 0, NX
        WRITE(10,'(2(1X,E13.6))') X(I), D(I)
      ENDDO
      WRITE(10,'(1X,E13.6,1X,I9)') DT,NT
      CLOSE(10)

      END


      SUBROUTINE STEP_FFT(D,XW,YW,DX,NX,NT,RX,TX,CT,RCT,TCT,DEN)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION CT(0:2*NX-1),RCT(0:2*NX-1),TCT(0:2*NX-1),DEN(0:2*NX-1)
      DIMENSION D(0:NX),XW(0:2*NX-1),YW(0:2*NX-1),XD(0:2*NX-1)
      DIMENSION YD(0:2*NX-1),XB(0:2*NX-1),YB(0:2*NX-1)
      DIMENSION XN(0:2*NX-1),YN(0:2*NX-1)

      DO I = 0, NX-1
        XD(I)=D(I); XD(NX+I)=D(NX-I); YD(I)=0.0; YD(NX+I)=0.0
      ENDDO
      DO I = 0, NX-1
        XB(I)=D(I)*D(I)*(-0.5+D(I)/3.); YB(I)=0.0
        XB(NX+I)=D(NX-I)*D(NX-I)*(-0.5+D(NX-I)/3.); YB(NX+I)=0.0
C        XB(I)=LOG(1.+D(I))-D(I); YB(I)=0.0
C        XB(NX+I)=LOG(1.+D(NX-I))-D(NX-I); YB(NX+I)=0.0
      ENDDO
      CALL MYFFT_FAST(XD,YD,2*NX,1)
      DO IT = 1, NT
        CALL MYFFT_FAST(XB,YB,2*NX,1)
        DO I = 0, 2*NX-1
          PX1=(XD(I)*XW(I)-YD(I)*YW(I))*RCT(I); PX2=XB(I)*TCT(I)+XD(I)
          PY1=(XD(I)*YW(I)+YD(I)*XW(I))*RCT(I); PY2=YB(I)*TCT(I)+YD(I)
          XD(I) = (PX1+PX2)/DEN(I); YD(I) = (PY1+PY2)/DEN(I)
          XN(I) = XD(I); YN(I) = YD(I)
        ENDDO
        CALL MYFFT_FAST(XN,YN,2*NX,-1)
        DO I = 0, 2*NX-1
          XB(I)=XN(I)*XN(I)*(-0.5+XN(I)/3.); YB(I)=0.0
C          XB(I)=LOG(1.+XN(I))-XN(I); YB(I)=0.0
        ENDDO
      ENDDO
      DO I = 0, NX
        D(I) = XN(I)
      ENDDO

      RETURN
      END


      SUBROUTINE STEP_FDM(D,W,XW,YW,DX,NX,DT,NT,RX,TX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(0:NX),W(0:NX),X(NX-1),B(NX-1),C(NX-1),E(NX-2),F(NX-2)
      DIMENSION VO(0:NX),XW(0:2*NX-1),YW(0:2*NX-1),FX(0:NX)

      RD = RX*DT/DX**2; TD = TX*DT/DX**2
      CI = -2.- 1./TD; RT = RX/TX
      DO I = 1, NX-2
        E(I) = 1; F(I) = 1
      ENDDO
      DO IT = 1, NT
        CALL CONVOL_FFT(D,XW,YW,VO,DX,NX)
C         CALL CONVOL_DIR(D,W,VO,DX,NX)
C     Neumann boundary conditions
        DO I = 0, NX
          FX(I) = D(I)*D(I)*(-0.5+D(I)/3.)
        ENDDO
        TEMA = 2.*(VO(1)-VO(0))
        TEMB = 2.*(FX(1)-FX(0))
        TEMC = 2.*(D(1)-D(0))
        D(0) = D(0) + RD*TEMA + TD*(TEMB+TEMC)
        TEMA = 2.*(VO(NX-1)-VO(NX))
        TEMB = 2.*(FX(NX-1)-FX(NX))
        TEMC = 2.*(D(NX-1)-D(NX))
        D(NX) = D(NX) + RD*TEMA + TD*(TEMB+TEMC)
        DO I = 1, NX-1
          C(I) = CI
          TED = -D(I)/TD; TEA = - RT*(VO(I+1)+VO(I-1)-VO(I)-VO(I))
          TEB = -(FX(I+1)+FX(I-1)-FX(I)-FX(I))
          IF(I.EQ.1) THEN
            B(I) = TED + TEA + TEB - D(0)
          ELSEIF(I.EQ.(NX-1)) THEN
            B(I) = TED + TEA + TEB - D(NX)
          ELSE
            B(I) = TED + TEA + TEB
          ENDIF
        ENDDO
        NG = NX - 1
        CALL CHASE(NG,C,E,F,X,B)
        DO I = 1, NX-1
          D(I) = X(I)
        ENDDO
      ENDDO

      RETURN
      END


      SUBROUTINE STEP_DIR(D,W,XW,YW,DX,NX,DT,NT,RX,TX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(0:NX),W(0:NX),XW(0:2*NX-1),YW(0:2*NX-1)
      DIMENSION VO(0:NX),FX(0:NX)

      RD = RX*DT/DX**2; TD = TX*DT/DX**2
      DO IT = 1, NT
c        CALL CONVOL_FFT(D,XW,YW,VO,DX,NX)
        CALL CONVOL_DIR(D,W,VO,DX,NX)
C     Neumann boundary conditions
        DO I = 0, NX
C          FX(I) = D(I)*(1.-0.5*D(I)+D(I)**3/3.)
          FX(I) = LOG(1.+D(I))
        ENDDO
        TEMA = 2.*(VO(1)-VO(0))
        TEMB = 2.*(FX(1)-FX(0))
        D(0) = D(0) + RD*TEMA + TD*TEMB
        TEMA = 2.*(VO(NX-1)-VO(NX))
        TEMB = 2.*(FX(NX-1)-FX(NX))
        D(NX) = D(NX) + RD*TEMA + TD*TEMB
        DO I = 1, NX-1
          TEMA = (VO(I+1)+VO(I-1)-2*VO(I))
          TEMB = (FX(I+1)+FX(I-1)-2*FX(I))
          D(I) =  D(I) + RD*TEMA + TD*TEMB
        ENDDO
        SUM = 0.0
        DO I = 0, NX
          IF((I.EQ.0).OR.(I.EQ.NX)) THEN
            SUM = SUM + 0.5*D(I)
          ELSE
            SUM = SUM + D(I)
          ENDIF
        ENDDO
      ENDDO

      RETURN
      END


      SUBROUTINE CHASE(N,C,D,E,X,B)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION X(N),B(N),C(N),D(N-1),E(N-1)
      DIMENSION F(N),G(N-1),U(N)

      F(1)=C(1)
      DO I=1,N-1
        G(I)=D(I)/F(I) 
        F(I+1)=C(I+1)-E(I)*G(I)
      ENDDO
      U(1)=B(1)/F(1)
      DO I=2,N
        U(I)=(B(I)-E(I-1)*U(I-1))/F(I) 
      ENDDO
      X(N)=U(N)
      DO I=N-1,1,-1
        X(I)=U(I)-G(I)*X(I+1)
      ENDDO

      RETURN
      END


      FUNCTION ENERGY(D,W,XW,YW,RX,TX,DX,NX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(0:NX),W(0:NX),VOL(0:NX),XW(0:2*NX-1),YW(0:2*NX-1)

C      CALL CONVOL_FFT(D,XW,YW,VOL,DX,NX)
      CALL CONVOL_DIR(D,W,VOL,DX,NX)
      SUM1 = 0.0; SUM2 = 0.0
      DO I = 0, NX
        IF((I.EQ.0).OR.(I.EQ.NX)) THEN
          SUM1 = SUM1 + 0.5*D(I)*VOL(I)
          SUM2 = SUM2 + 0.5*(1.+D(I))*LOG(1.+D(I))
        ELSE
          SUM1 = SUM1 + D(I)*VOL(I)
          SUM2 = SUM2 + (1.+D(I))*LOG(1.+D(I))
        ENDIF
      ENDDO
      SUM1 = SUM1*RX**2/2.; SUM2 = SUM2*RX*TX
      ENERGY = (SUM1 + SUM2)/NX

      RETURN
      END

      SUBROUTINE CONVOL_FFT(D,XW,YW,VOL,DX,NX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (PI=3.14159265358979323845)
      DIMENSION D(0:NX),W(0:NX),VOL(0:NX),XD(0:2*NX-1),YD(0:2*NX-1)
      DIMENSION XW(0:2*NX-1),YW(0:2*NX-1),XV(0:2*NX-1),YV(0:2*NX-1)

      DO I = 0, NX-1
        XD(I) = D(I); XD(NX+I) = D(NX-I); YD(I) = 0.0; YD(NX+I) = 0.0
      ENDDO
      CALL MYFFT_FAST(XD,YD,2*NX,1)
      DO I = 0, 2*NX-1
        XV(I) = XD(I)*XW(I) - YD(I)*YW(I)
        YV(I) = XD(I)*YW(I) + XW(I)*YD(I)
      ENDDO
      CALL MYFFT_FAST(XV,YV,2*NX,-1)     
      DO I = 0, NX
        VOL(I) = XV(I)*DX
      ENDDO

      RETURN
      END


      SUBROUTINE CONVOL_DIR(D,W,VOL,DX,NX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(0:NX),W(0:NX),VOL(0:NX)

      DO I = 0, NX
        SUM = 0.0
        DO J = -NX, NX
          NJ = ABS(J); K = ABS(I+J)
          IF(K.GT.NX) K = 2*NX-K
          IF(NJ.EQ.NX) THEN
            SUM = SUM + 0.5*W(NJ)*D(K)
          ELSE
            SUM = SUM + W(NJ)*D(K)
          ENDIF
        ENDDO
        VOL(I) = SUM*DX
      ENDDO

      RETURN
      END

      SUBROUTINE FFTVAL(CT,RCT,TCT,DEN,DX,NX,DT,RX,TX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(PI=3.1415926535897932384626)
      DIMENSION CT(0:2*NX-1),RCT(0:2*NX-1),TCT(0:2*NX-1),DEN(0:2*NX-1)
      
      DC = PI/(NX*DX)
      DO I = 0, NX-1
        CS1 = I*DC; CS2 = (I-NX)*DC
        CT(I) = -DT*CS1**2; CT(NX+I)=-DT*CS2**2
        RCT(I) = RX*CT(I)*DX; RCT(NX+I) = RX*CT(NX+I)*DX
        TCT(I) = TX*CT(I); TCT(NX+I) = TX*CT(NX+I)
        DEN(I) = 1.-TCT(I); DEN(NX+I) = 1.-TCT(NX+I)
      ENDDO

      RETURN
      END

      SUBROUTINE WVAL(W,WX,WY,XIN,XOUT,DX,NX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER(PI=3.1415926535897932384626)
      DIMENSION W(0:NX),WX(0:2*NX-1),WY(0:2*NX-1)
      
      DO I = 0, NX
        XI = I*DX; W(I) = WLJ(XI,XIN,XOUT)
      ENDDO
      DO I = 0, NX-1
        WX(I) = W(I); WX(NX+I) = W(NX-I)
        WY(I) = 0.0; WY(NX+I) = 0.0
      ENDDO
      CALL MYFFT_FAST(WX,WY,2*NX,1)

      RETURN
      END


      FUNCTION WLJ(X,XIN,XOUT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      
      WLJC = 1./(XOUT**12) - 2./(XOUT**6)
      A = -6.*(1./XIN**6-1.)/XIN**8; B = (7./XIN**6-8.)/XIN**6
      IF(X.LE.XIN) THEN
        WLJ =  A*X**2+B - WLJC
      ELSEIF(X.GE.XOUT) THEN
        WLJ = 0.0
      ELSE
        WLJ =  1./(X**12) - 2./(X**6) - WLJC
      ENDIF

      RETURN
      END


      SUBROUTINE MYFFT_FAST(X,Y,N,IDR)
       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION X(*), Y(*)
      PARAMETER(TWOPI = 6.2831853071795864769)

      IF ( N .EQ. 1 ) RETURN
      IF ( IDR .EQ. -1 ) THEN
        DO I = 1, N
          Y(I) = - Y(I)
        ENDDO
      ENDIF
      N2 = 2 * N
 
      E = TWOPI / N2
      CE1 = 2.*SIN(0.5*E)**2; CI1 = SIN(E)
      CE3 = 2.*SIN(1.5*E)**2; CI3 = SIN(3.*E)
      DO WHILE(N2.GE.4)
        N2 = N2 / 2;  N4 = N2 / 4
        TE1 = CE1; TE3 = CE3
        CE1 = 2.*CI1**2; CI1 = 2.*(CI1-CI1*TE1)
        CE3 = 2.*CI3**2; CI3 = 2.*(CI3-CI3*TE3)
        CC1 = 1.0; CC3 = 1.0; SS1 = 0.0; SS3 = 0.0
        DO J = 1, N4
          IS = J; ID = 2 * N2
          DO WHILE(IS.LT.N)
            DO I0 = IS, N-1, ID
              I1 = I0 + N4; I2 = I1 + N4; I3 = I2 + N4
              R1 = X(I0) - X(I2); X(I0) = X(I0) + X(I2)
              R2 = X(I1) - X(I3); X(I1) = X(I1) + X(I3)
              S1 = Y(I0) - Y(I2); Y(I0) = Y(I0) + Y(I2)
              S2 = Y(I1) - Y(I3); Y(I1) = Y(I1) + Y(I3)
              S3 = R1 - S2; R1 = R1 + S2
              S2 = R2 - S1; R2 = R2 + S1
              X(I2) = R1 * CC1 - S2 * SS1
              Y(I2) = - S2 * CC1 - R1 * SS1
              X(I3) = S3 * CC3 + R2 * SS3
              Y(I3) = R2 * CC3 - S3 * SS3
            ENDDO
            IS = 2 * ID - N2 + J
            ID = 4 * ID
          ENDDO
          TC1 = CC1; TC3 = CC3
          CC1 = CC1 - CC1*CE1 - SS1*CI1; SS1 = SS1 - SS1*CE1 + TC1*CI1
          CC3 = CC3 - CC3*CE3 - SS3*CI3; SS3 = SS3 - SS3*CE3 + TC3*CI3
        ENDDO
      ENDDO
 
      IS = 1; ID = 4
      DO WHILE(IS.LT.N)
        DO I0 = IS, N, ID
          I1 = I0 + 1;
          R1 = X(I0); X(I0) = R1 + X(I1); X(I1) = R1 - X(I1)
          R1 = Y(I0); Y(I0) = R1 + Y(I1); Y(I1) = R1 - Y(I1)
        ENDDO
        IS = 2 * ID - 1
        ID = 4 * ID
      ENDDO
 
      J = 1; N1 = N - 1
      DO I = 1, N1
        IF(I.LT.J) THEN
          XT = X(J); X(J) = X(I)
          X(I) = XT; XT = Y(J)
          Y(J) = Y(I); Y(I) = XT
        ENDIF
        K = N / 2
        DO WHILE(K.LT.J)
          J = J -K; K = K/2
        ENDDO
        J = J + K
      ENDDO

      IF (IDR.EQ. -1 ) THEN
        DO I = 1, N
          X(I) = X(I) / N; Y(I) = - Y(I) / N
        ENDDO
      ENDIF

      RETURN
      END
