      !-------------------Programma Diplomatikhs Ergasias - Part 2 - ----------------
  

      MODULE SINUS_SOLVER
      
      IMPLICIT NONE     
                                                          !P0  = Extra chordwise C.P
      INTEGER :: NR_MAX ,NCH_MAX ,MOMENTS_MAX ,M0,N0 ,P0  !CPW = number of spanwise C.P on wake C.P 
      INTEGER :: MNT_MAX ,NBL_MAX ,MNT ,CPW ,CPB          !CPB = number of control points on blade
      PARAMETER(NR_MAX=13,NCH_MAX=17,NBL_MAX=1,MOMENTS_MAX=200,M0=9,N0=7,P0=0)   
      PARAMETER(MNT_MAX=M0*(N0+1)*NBL_MAX, CPW=10*(M0+2-1)+1 ) 
      PARAMETER(CPB=M0*(N0+1+P0)*NBL_MAX) 
      INTEGER :: I,J,K,L,M,N,Z ,II,JJ,KZ,KZZ,IN  ,M1,N1,MID    
      INTEGER :: IE0,IE1,IE2,IOPT,IREG ,LL(MNT_MAX) 
 
     !Variables for INPUT_GEOMETRY
      INTEGER :: NBL,NR,NCH,MOMENTS
      DOUBLE PRECISION :: ROUT ,RHUB
      DOUBLE PRECISION,DIMENSION(NR_MAX) :: U11 ,XS1,DXS1 ,SKEW1,DSKEW1 ,FI1,DFI1 ,CHORD1,DCHORD1 
      DOUBLE PRECISION,DIMENSION(NR_MAX,NCH_MAX):: U21 ,DEC_DU1,DEC_DU2  ,DET_DU1,DET_DU2
      DOUBLE PRECISION :: T(MOMENTS_MAX),DT
      CHARACTER(LEN=21) :: MY_LABEL
     !Variables for INPUT_POSITION_VELOCITY
      DOUBLE PRECISION,DIMENSION(NR_MAX,NCH_MAX,NBL_MAX,3) :: FACE  ,BACK ,VFACE  ,VBACK
      DOUBLE PRECISION ,DIMENSION(3) :: ORIGIN,UNIT1,UNIT2,UNIT3
     !Variables for WAKE_POSITION, WAKE_INERTIAL, WAKE_DIPOLE_INTENSITY ,WAKE_CONTRIBUTION 
      DOUBLE PRECISION :: WAKEP(CPW,NBL_MAX,MOMENTS_MAX,3) 
      DOUBLE PRECISION :: WAKEP_I(CPW,NBL_MAX,MOMENTS_MAX,3),DR(3)
      DOUBLE PRECISION :: M_WAKE(CPW,NBL_MAX,MOMENTS_MAX),WP(3)  
     !Variables for INDUCTION_FACTOR
      DOUBLE PRECISION ,DIMENSION(M0,3)    :: BSM0 ,CSM0  ,SS_S0 ,SS_S01 ,SS_S00 ,SS_S02 
      DOUBLE PRECISION ,DIMENSION(M0,N0,3) :: BSSMN,CSSMN ,SS_SSN,SS_SSN1,SS_SSN0,SS_SSN2
      DOUBLE PRECISION ,DIMENSION(M0+2,P0+N0+1,NBL_MAX,M0,3)    :: TS_M0    
      DOUBLE PRECISION ,DIMENSION(M0+2,P0+N0+1,NBL_MAX,M0,N0,3) :: TSS_MN
      DOUBLE PRECISION :: THT(M0+2,P0+N0+1,NBL_MAX,MNT_MAX,3) ,THTP(MNT_MAX,3)
      DOUBLE PRECISION ,DIMENSION(3) :: E,EC0
      DOUBLE PRECISION :: EBS(M0,3),ECS(N0,3)
      DOUBLE PRECISION :: F0,F1,AB1,AB2,Y1A,Y2A
     !Variables for LINEAR_SYSTEM_FORMATION
      DOUBLE PRECISION :: VT(CPB,MNT_MAX),B(CPB) ,A(MNT_MAX,MOMENTS_MAX+1) 
      DOUBLE PRECISION :: P, Q_SM0,Q_SSMN 
     !Variables for  PRESSURE_DISTRIBUTION, LIFT_THRUST    
      DOUBLE PRECISION :: AA11(P0+N0+2),AA12(M0+2),DPRE(M0+2,P0+N0+2,NBL_MAX)
      DOUBLE PRECISION ,DIMENSION(M0+2,NBL_MAX) :: GAMMA, LE_SF, DL, Q3_PHYSICAL
      DOUBLE PRECISION :: THRUST(MOMENTS_MAX,3),THRUST_I(MOMENTS_MAX,3),LIFT(NBL_MAX)   
     
     !Variables for PROPERTIES_AT_CP
     !Auxilary variables for the interpolation of input data  at C.P
      DOUBLE PRECISION ,DIMENSION(NR_MAX) :: AA1,AA2,AA3,AA4,AA5,AA6,AA7,AA8,AA9,AA10
      DOUBLE PRECISION ,DIMENSION(NCH_MAX):: BB1,BB2,BB3,BB4,BB5,BB6,BB7,BB8,BB9,BB10,BB11       
     !Variables for Geometric Properties at C.P
      DOUBLE PRECISION :: YP1,YP2,UP1,UP2 ,CHORDP ,XSP,SKEWP,FIP,DCHORDP,DXSP,DSKEWP,DFIP
      DOUBLE PRECISION :: DELTAP,XP1,XP2,XP3,ARG 
      DOUBLE PRECISION :: G11P,G21P,G31P,G12P,G22P,G32P,G13P,G23P,G33P ,G 
      DOUBLE PRECISION :: C11P,C12P,C13P,C21P,C22P,C23P,C31P,C32P,C33P ,GP
      DOUBLE PRECISION :: GG13,GG33,GG22 !GG,GG11,GG21,GG12,GG32,GG31,GG23,C
      DOUBLE PRECISION :: A11P,A22P,A21P ,ALPHAP
      DOUBLE PRECISION :: DEC_DU1P,DEC_DU2P,DET_DU1P,DET_DU2P, NP(M0+2,P0+N0+2,3)
      DOUBLE PRECISION :: E2(M0+2,NBL_MAX,3)
     !Variables for Kinematic Properties at C.P 
      DOUBLE PRECISION ,DIMENSION(3) :: VFACEP,VBACKP,WFACEP,WBACKP
      DOUBLE PRECISION :: QQ(M0+2,P0+N0+2,NBL_MAX,4),Q(4)
      DOUBLE PRECISION ,DIMENSION(M0+2,P0+N0+2,NBL_MAX):: S ,DPRE_S
     !Variables for Geometric Properties at Integration Point
      DOUBLE PRECISION :: Y1,Y2,U1,U2,CHR,XS,SK,FI,DCHR,DXS,DSK,DFI
      DOUBLE PRECISION :: DELTA,X1,X2,X3
      DOUBLE PRECISION :: G11,G21,G31,G12,G22,G32,G13,G23,G33
      
      DOUBLE PRECISION :: SPAN ,PI ,RII ,RJJ ,EE ,D
      !Test Variables
      DOUBLE PRECISION :: VT1(CPB,MNT_MAX),B1(CPB)

      
      CONTAINS


!     ******************************************************************
!     *                                                                *
!     *                        SUBROUTINE MAIN                         *
!     *                                                                *
!     *                FOR THE CALCULATION OF A_{MNZ,IN}               *
!     *                                                                *
!     ******************************************************************

      SUBROUTINE MAIN

    ! introduce pi
      PI = 4.D0*ATAN(1.D0)       
   
      OPEN (unit=3 ,file='Input_for_Solver.txt')
      OPEN (unit=4 ,file='Propeller_Graph.txt')
      OPEN (unit=5 ,file='My_output_Unstable.txt')
      
      CALL INPUT_GEOMETRY

      M1 = M0+1           !M1+1 C.P in the y1 direction including those on tip and hub   
      N1 = P0+N0+1        !N1+1 C.P in the y2 direction including those on L.E and T.E                 
      SPAN = ROUT-RHUB    !Equal spacing of C.P in the Y1,Y2 coordinate system 
      MNT  = M0*(N0+1)*NBL
       
     
      DO IN=1,MOMENTS !IN=Instant 
     !The first instnace (IN=1) the blade has full velocity and is at its starting position (step velocity) 
     !The instant before that (IN=0) the blade has zero velocity and is at its starting position      
                
        !Initialization of system is crucial since therea are some UNIDENTIFIED 
        !cells in the Bernolli equation (e.g. M=2,N=15 ->K=99 defines L=1,15,...)
        !DGAUSS and the pressure calculation changes the initial
        DO I=1,CPB ;DO J=1,MNT_MAX ;VT(I,J) = 0.D0;END DO ;END DO 
        DO I=1,CPB ;B(I) = 0.D0;END DO

        CALL INPUT_POSITION_VELOCITY  
        CALL WAKE_POSITION
        CALL WAKE_INERTIAL
        CALL WAKE_DIPOLE_INTENSITY

       ! ----------APPROXIMATION OF VORTICITY WITH DOUBLE SIN-SERIES.----------------
        DO Z=1,NBL 
        DO M=2,M1 
           !For M=1-->YP1=0 (hub) and M=M1+1-->YP1=PI (tip)
           !For the no-entrance b.c you need 0<Y1P<PI --> M=2,M1
            YP1 = PI/DBLE(M1)*DBLE(M-1) 
            UP1 = 0.5D0*(ROUT+RHUB) - 0.5D0*(ROUT-RHUB)*COS(YP1)
           !UP1 = SPAN/DBLE(M1)*DBLE(M-1)
           !IF (UP1>SPAN) THEN; UP1=SPAN;END IF
           !YP1 = ACOS((0.5D0*(ROUT+RHUB)-UP1)/(0.5D0*(ROUT-RHUB)))

            CHORDP  = FILLIN(UP1,U11,CHORD1,NR) 
            XSP     = FILLIN(UP1,U11,XS1,NR)
            SKEWP   = FILLIN(UP1,U11,SKEW1,NR)
            FIP     = FILLIN(UP1,U11,FI1,NR)
            DCHORDP = FILLIN(UP1,U11,DCHORD1,NR)
            DXSP    = FILLIN(UP1,U11,DXS1,NR)
            DSKEWP  = FILLIN(UP1,U11,DSKEW1,NR)
            DFIP    = FILLIN(UP1,U11,DFI1,NR)
        DO N=1,N1+1 !Point N1+1 refers to the L.E
            WRITE(*,25) M,N,Z,IN
         25 FORMAT(' COMPUTATION AT CONTROL POINT P(',I2,',',I2,',',I2,') ,INSTANT =',I4)  
           !For N=0(N1+1)-->YP2=0 (L.E) and N=N1-->YP1=PI (T.E)
           !For the no-entrance b.c you need 0<Y2P<PI --> N=1,N1-1
            IF (N<N1+1) THEN
            YP2   = PI/DBLE(N1)*DBLE(N)
            ELSE !N==N1+1
            YP2   = 0.D0
            END IF     
            UP2   =-0.5D0*COS(YP2)
            !UP2   = 0.5*(0.3-0.5) - 0.5*(0.3+0.5)*COS(YP2)
            CALL PROPERTIES_AT_CP            
            CALL LINEAR_SYSTEM_FORMATION
        END DO !DO N=1,N1+1          
        END DO !DO M=1,M1+1 
        END DO !DO Z=1,NBL 
        
        DO I=1,CPB ;DO J=1,MNT_MAX ;VT1(I,J) = 0.D0 ;END DO ;END DO 
        DO I=1,CPB ;DO J=1,MNT ;VT1(I,J) = VT(I,J) ;END DO ;END DO          
        DO I=1,CPB ; B1(I) = B(I) ;END DO
         
        !IF (IN==35) THEN   
        !DO I=1,MNT 
        !  WRITE(5,*) '-----------------I=',I,'-------------'
        !  DO J=1,MNT
        !  WRITE(5,*) 'VT1(',I,',',J,')=', VT1(I,J)
        !END DO ;END DO
        !END IF
        !DO I=1,MNT 
        !  WRITE(5,'(56F6.1)') (VT1(I,J),J=1,K)
        !END DO 
        !DO I=1,MNT
        !  WRITE(5,*) B1(I),'---'
        !END DO 
        !IF (IN==5) THEN
        !close(5)
        !END IF
        
        IOPT=1
        CALL DGAUSS(MNT_MAX,MNT,VT,B,IOPT,D,IREG,LL)
        IOPT=2
        CALL DGAUSS(MNT_MAX,MNT,VT,B,IOPT,D,IREG,LL) 
        ! From this point on , vector B is the solution of the system

        !WRITE(5,*) '-------INSTANT=',IN,'------------'
        DO L=1,MNT               
          A(L,IN)= B(L)
          !WRITE(5,*) 'B(',L,')=',B(L) 
          !WRITE(5,*) B(L) 
        END DO    
        !WRITE(5,*) IN ,B(2)      
        DO I=1,MNT 
         B(I) = 0.D0
         DO J=1,MNT
           !WRITE(5,*) 'VT1(',I,',',J,')=', VT1(I,J)
           B(I) = B(I)+A(J,IN)*VT1(I,J)
         END DO
         !IF (ABS(B(I)-B1(I))>0.00001D0) THEN
            WRITE(*,*) B(I)-B1(I) ,'yo'
         !END IF
        END DO

        CALL PRESSURE_DISTRIBUTION
        CALL LIFT_THRUST
                           
      END DO !DO IN=1,MOMENTS
       
      END SUBROUTINE MAIN


!     ******************************************************************
!     *                                                                *
!     *                  SUBROUTINE INPUT_GEOMETRY                     *
!     *                                                                *
!     *            INPUT OF BASIC GEOMETRIC AND TIME STEPS             *
!     *                                                                *
!     ****************************************************************** 
  
      SUBROUTINE INPUT_GEOMETRY     
      INTEGER :: I
                          
      READ(3,*)  ROUT ,RHUB ,NBL ,NR ,NCH ,MOMENTS 
      WRITE(4,6) M0 ,N0 ,NBL ,NR ,NCH ,MOMENTS 
    6 FORMAT (6I3)
      DO I=1,NR
      READ(3,*) U11(I), CHORD1(I),XS1(I),SKEW1(I),FI1(I), DCHORD1(I),DXS1(I),DSKEW1(I),DFI1(I)
      END DO
      DO I=1,NR
        READ(3,15) (U21(I,J),J=1,NCH),(DEC_DU1(I,J),J=1,NCH),(DEC_DU2(I,J),J=1,NCH)
        READ(3,15) (DET_DU1(I,J),J=1,NCH) ,(DET_DU2(I,J),J=1,NCH)
      END DO
   15 FORMAT (11F7.4/6F7.4)

      READ(3,'(A7)') MY_LABEL !'Moments'
      DO I=1,MOMENTS
        READ(3,'(F7.4)') T(I) !For constant timestep we need only the MOMENTS and DT 
      END DO
      DT = T(2)-T(1)    !We assume constant timesteps   

      END SUBROUTINE INPUT_GEOMETRY


      
!     ******************************************************************
!     *                                                                *
!     *               SUBROUTINE INPUT POSITION VELOCITY               *
!     *                                                                *
!     *          INPUT OF POSITION,VELOCITIES AND UNIT VECTORS         *
!     *                                                                *
!     ******************************************************************

      SUBROUTINE INPUT_POSITION_VELOCITY

      INTEGER :: I,J,K

        READ(3,'(A20)') MY_LABEL !'Position Vector FACE'
        DO K=1,NBL ;DO I=1,NR ;DO J=1,NCH
          READ(3,20) FACE(I,J,K,1) ,FACE(I,J,K,2) ,FACE(I,J,K,3)
          WRITE(4,20) FACE(I,J,K,1) ,FACE(I,J,K,2) ,FACE(I,J,K,3)
        END DO ;END DO ;END DO
        READ(3,'(A20)') MY_LABEL !'Position Vector BACK'
        DO K=1,NBL ;DO I=1,NR ;DO J=1,NCH
          READ(3,20) BACK(I,J,K,1) ,BACK(I,J,K,2) ,BACK(I,J,K,3)
          WRITE(4,20) BACK(I,J,K,1) ,BACK(I,J,K,2) ,BACK(I,J,K,3)
        END DO ;END DO ;END DO   
        READ(3,'(A20)') MY_LABEL !'Position Vector Path' 
        READ(3,20) ORIGIN(1),ORIGIN(2),ORIGIN(3)
        READ(3,'(A20)') MY_LABEL !'Unit Vector 1'
        READ(3,20) UNIT1(1),UNIT1(2),UNIT1(3)
        READ(3,'(A20)') MY_LABEL !'Unit Vector 1'
        READ(3,20) UNIT2(1),UNIT2(2),UNIT2(3)
        READ(3,'(A20)') MY_LABEL !'Unit Vector 1'
        READ(3,20) UNIT3(1),UNIT3(2),UNIT3(3)
     20 FORMAT(3F16.5)
        READ(3,'(A15)') MY_LABEL !'Velocities face'  
        DO K=1,NBL ;DO I=1,NR ;DO J=1,NCH
          READ(3,20) VFACE(I,J,K,1) ,VFACE(I,J,K,2) ,VFACE(I,J,K,3)
        END DO ;END DO ;END DO
        READ(3,'(A15)') MY_LABEL !'Velocities back'
        DO K=1,NBL ;DO I=1,NR ;DO J=1,NCH
          READ(3,20) VBACK(I,J,K,1) ,VBACK(I,J,K,2) ,VBACK(I,J,K,3)
        END DO ;END DO ;END DO
        
      END SUBROUTINE INPUT_POSITION_VELOCITY


      
!     ******************************************************************
!     *                                                                *
!     *                     SUBROUTINE WAKE POSITION                   *
!     *                                                                *
!     *  ADDS THE KUTTA STRIP AND EVERY WAKE STRIP ONE AT A TIME STEP  *
!     *                                                                *
!     ******************************************************************
 
      SUBROUTINE WAKE_POSITION
      
      INTEGER :: Z,M      

       ! Calculate the position vectors of every point on the wake 
       ! so far ,using the inertial coordinate system
        DO Z=1,NBL 
        DO M=1,CPW   
         !For M=1-->YP1=0 (hub) and M=CPW-->YP1=PI (tip)
         !Calculation of the C.P of the wake

         IF ((IN==1).OR.(IN==2)) THEN
         !Notice that the wake panels smaller near the tip and hub because equal distributed 
         !C.P in the YP1 direction are not equally distributed in the UP1 direction.But that
         !is ok because the vortex intesifies near the tip and hub so we need finer discretization
          YP1 = PI/DBLE(CPW-1)*DBLE(M-1) 
          UP1 = 0.5D0*(ROUT+RHUB) - 0.5D0*(ROUT-RHUB)*COS(YP1)
          DELTAP = 2.D0*PI*(DBLE(Z)-1.D0)/DBLE(NBL)
            
          CHORDP  = FILLIN(UP1,U11,CHORD1,NR)
          XSP     = FILLIN(UP1,U11,XS1,NR)
          SKEWP   = FILLIN(UP1,U11,SKEW1,NR)
          FIP     = FILLIN(UP1,U11,FI1,NR)
          DCHORDP = FILLIN(UP1,U11,DCHORD1,NR)
          DXSP    = FILLIN(UP1,U11,DXS1,NR)
          DSKEWP  = FILLIN(UP1,U11,DSKEW1,NR)
          DFIP    = FILLIN(UP1,U11,DFI1,NR)
          IF (IN==1) THEN
          UP2 = 0.5D0 
          CALL GEOM(NBL,UP1,UP2,DELTAP,CHORDP,XSP,SKEWP,FIP,DCHORDP,DXSP,DSKEWP,DFIP,& 
                    XP1,XP2,XP3,ARG,G11P,G21P,G31P,G12P,G22P,G32P,G13P,G23P,G33P)
         !Position Vector of C.P on wake ,projected on the moving cartesian system           
          WAKEP(M,Z,IN,1)=XP1 ;WAKEP(M,Z,IN,2)=XP2 ;WAKEP(M,Z,IN,3)=XP3

          ELSE !(IN==2)
         !Position Vector of C.P on wake ,projected on the moving cartesian system
          WAKEP(M,Z,IN,1)=WAKEP(M,Z,IN-1,1) 
          WAKEP(M,Z,IN,2)=WAKEP(M,Z,IN-1,2) 
          WAKEP(M,Z,IN,3)=WAKEP(M,Z,IN-1,3)          
          UP2 = 0.5D0 + 0.1D0
          CALL GEOM(NBL,UP1,UP2,DELTAP,CHORDP,XSP,SKEWP,FIP,DCHORDP,DXSP,DSKEWP,DFIP,& 
                    XP1,XP2,XP3,ARG,G11P,G21P,G31P,G12P,G22P,G32P,G13P,G23P,G33P)
         !Position Vector of C.P on wake ,projected on the moving cartesian system
          WAKEP(M,Z,IN-1,1)=XP1 ;WAKEP(M,Z,IN-1,2)=XP2 ;WAKEP(M,Z,IN-1,3)=XP3
          END IF !IF (IN==1)
          
          ELSE !(IN>=3)
         !Position Vector of C.P on wake ,projected on the moving cartesian system
          WAKEP(M,Z,IN,1)=WAKEP(M,Z,IN-1,1) 
          WAKEP(M,Z,IN,2)=WAKEP(M,Z,IN-1,2) 
          WAKEP(M,Z,IN,3)=WAKEP(M,Z,IN-1,3)

         !Position Vector of C.P on wake ,projected on the moving cartesian system
          WAKEP(M,Z,IN-1,1)=WAKEP(M,Z,IN-2,1) 
          WAKEP(M,Z,IN-1,2)=WAKEP(M,Z,IN-2,2) 
          WAKEP(M,Z,IN-1,3)=WAKEP(M,Z,IN-2,3)
          END IF
        END DO !DO M=1,CPW
        END DO !DO Z=1,NBL 

        END SUBROUTINE WAKE_POSITION      



!     ******************************************************************
!     *                                                                *
!     *                    SUBROUTINE WAKE_INERTIAL                    *
!     *                                                                *
!     *    CALCULATES THE POSITION OF WAKE IN THE INERTIAL SYSTEM      *
!     *                                                                *
!     ******************************************************************

      SUBROUTINE WAKE_INERTIAL

      INTEGER :: I,Z,M
       
        DO I=1,IN
        DO Z=1,NBL
        DO M=1,CPW
          IF ((I==IN-1).OR.(I==IN)) THEN
          !Position of T.E and Kutta strip in the inertial system expressed with the inertial system's components  
           WAKEP_I(M,Z,I,1) =WAKEP(M,Z,I,1)*UNIT1(1) +WAKEP(M,Z,I,2)*UNIT2(1) +WAKEP(M,Z,I,3)*UNIT3(1)   !----------------------------------------------------------------------------------Put these in a do-loop (DO J=1,3).You have to make an array out of UNIT1,UNIT2,UNIT3
           WAKEP_I(M,Z,I,2) =WAKEP(M,Z,I,1)*UNIT1(2) +WAKEP(M,Z,I,2)*UNIT2(2) +WAKEP(M,Z,I,3)*UNIT3(2)
           WAKEP_I(M,Z,I,3) =WAKEP(M,Z,I,1)*UNIT1(3) +WAKEP(M,Z,I,2)*UNIT2(3) +WAKEP(M,Z,I,3)*UNIT3(3)  

           WAKEP_I(M,Z,I,1) =WAKEP_I(M,Z,I,1) + ORIGIN(1)
           WAKEP_I(M,Z,I,2) =WAKEP_I(M,Z,I,2) + ORIGIN(2)
           WAKEP_I(M,Z,I,3) =WAKEP_I(M,Z,I,3) + ORIGIN(3)
          ELSE !I<=IN-2
          !Position of wake in the moving cartesian system expressed with the moving system's components
           WAKEP(M,Z,I,1)=WAKEP_I(M,Z,I,1)*UNIT1(1)+WAKEP_I(M,Z,I,2)*UNIT1(2)+WAKEP_I(M,Z,I,3)*UNIT1(3) 
           WAKEP(M,Z,I,2)=WAKEP_I(M,Z,I,1)*UNIT2(1)+WAKEP_I(M,Z,I,2)*UNIT2(2)+WAKEP_I(M,Z,I,3)*UNIT2(3)
           WAKEP(M,Z,I,3)=WAKEP_I(M,Z,I,1)*UNIT3(1)+WAKEP_I(M,Z,I,2)*UNIT3(2)+WAKEP_I(M,Z,I,3)*UNIT3(3) 
          !Parallel transition of origin point projected on the moving cartesian system
           DR(1) = ORIGIN(1)*UNIT1(1) + ORIGIN(2)*UNIT1(2) + ORIGIN(3)*UNIT1(3)
           DR(2) = ORIGIN(1)*UNIT2(1) + ORIGIN(2)*UNIT2(2) + ORIGIN(3)*UNIT2(3)
           DR(3) = ORIGIN(1)*UNIT3(1) + ORIGIN(2)*UNIT3(2) + ORIGIN(3)*UNIT3(3)
          !Relative (minus sign) parallel transition of wake
           WAKEP(M,Z,I,1) = WAKEP(M,Z,I,1) - DR(1)
           WAKEP(M,Z,I,2) = WAKEP(M,Z,I,2) - DR(2)
           WAKEP(M,Z,I,3) = WAKEP(M,Z,I,3) - DR(3) 
          END IF  
         !Wake as seen by an inertial observer
          WRITE(4,20) WAKEP_I(M,Z,I,1) ,WAKEP_I(M,Z,I,2) ,WAKEP_I(M,Z,I,3)
       20 FORMAT(3F16.5)
         !Wake as seen by an observer on the blades
          !WRITE(4,20) WAKEP(M,Z,I,1) ,WAKEP(M,Z,I,2) ,WAKEP(M,Z,I,3)  
        END DO !DO M=1,CPW
        END DO !DO Z=1,NBL 
        END DO !DO I=1,IN  
            
      END SUBROUTINE WAKE_INERTIAL      


!     ******************************************************************
!     *                                                                *
!     *               SUBROUTINE WAKE_DIPOLE_INTENSITY                 *
!     *                                                                *
!     *    CALCULATES THE DIPOLE INTENSITY AT EACH NEW WAKE STRIP      *
!     *                                                                *
!     ******************************************************************

      SUBROUTINE WAKE_DIPOLE_INTENSITY

      INTEGER :: Z,M,II,JJ,L

        !Calculation of dipole intesity on the C.P of wake
        IF (IN>=2) THEN 
        DO Z=1,NBL 
        DO M=1,CPW
            YP1 = PI/DBLE(CPW-1)*DBLE(M-1) 
            M_WAKE(M,Z,IN-1) = 0.D0

           !For the Sinus Series, mu=0 at tip and hub            
            IF ((M>1).AND.(M<CPW)) THEN             
            DO II=1,M0
              RII = DBLE(II)         
              L = II + (Z-1)*M0*(N0+1)
              M_WAKE(M,Z,IN-1) = M_WAKE(M,Z,IN-1) + A(L,IN-1)*SIN(RII*YP1)*PI/2.D0 !A_SM0
              DO JJ=1,N0
                RJJ = DBLE(JJ)
                L = II + M0*JJ + (Z-1)*M0*(N0+1)
                M_WAKE(M,Z,IN-1) = M_WAKE(M,Z,IN-1) + A(L,IN-1)*SIN(RII*YP1)*(1.D0-COS(RJJ*PI))/RJJ !A_SSMN                 
            END DO ;END DO            
            END IF          
        END DO !DO M=1,CPW
        END DO !DO Z=1,NBL
        END IF       

      END SUBROUTINE WAKE_DIPOLE_INTENSITY      

      

!     ******************************************************************
!     *                                                                *
!     *                   SUBROUTINE PROPERTIES_AT_CP                  *
!     *                                                                *
!     *    CALCULATES THE GEOMETRIC AND KINEMATIC PROPERTIES AT C.P    *
!     *                                                                *
!     ******************************************************************

      SUBROUTINE PROPERTIES_AT_CP

      INTEGER :: I,J 

      DELTAP= 2.D0*PI*(DBLE(Z)-1.D0)/DBLE(NBL)  !The C.P belong on the reference blade    
      CALL GEOM(NBL,UP1,UP2,DELTAP,CHORDP,XSP,SKEWP,FIP,DCHORDP,DXSP,DSKEWP,DFIP,& 
      XP1,XP2,XP3,ARG,G11P,G21P,G31P,G12P,G22P,G32P,G13P,G23P,G33P)  
     !Calculation of surface metric alpha
      A11P = G11P**2 + G21P**2 + G31P**2
      A22P = G12P**2 + G22P**2 + G32P**2
      A21P = G12P*G11P + G22P*G21P + G32P*G31P
      ALPHAP = A11P*A22P - A21P**2 
            
     !Just to be sure we calculate the metric G at CP analyticaly
     !GG11 = G11P*G11P + G21P*G21P + G31P*G31P
     !GG12 = G11P*G12P + G21P*G22P + G31P*G32P
      GG13 = G11P*G13P + G21P*G23P + G31P*G33P
     !GG21 = G12P*G11P + G22P*G21P + G32P*G31P
      GG22 = G12P*G12P + G22P*G22P + G32P*G32P
     !GG23 = G12P*G13P + G22P*G23P + G32P*G33P
     !GG31 = G13P*G11P + G23P*G21P + G33P*G31P
     !GG32 = G13P*G12P + G23P*G22P + G33P*G32P
      GG33 = G13P*G13P + G23P*G23P + G33P*G33P
     !GP=GG11*(GG22*GG33-GG23*GG32)-GG12*(GG21*GG33-GG23*GG31)+GG13*(GG21*GG32-GG22*GG31)
     !WRITE(*,*) 'GP=',GP,'CHORD^4=',CHORDP**4 
      GP = CHORDP**4 

     !Interpolation of input data for C.P
      DO I=1,NR 
        DO J=1,NCH 
          BB1(J) = U21(I,J)    
          BB2(J) = DEC_DU1(I,J) 
          BB3(J) = DEC_DU2(I,J)
          BB4(J) = DET_DU1(I,J)
          BB5(J) = DET_DU2(I,J)
          BB6(J) = VFACE(I,J,Z,1)
          BB7(J) = VFACE(I,J,Z,2)
          BB8(J) = VFACE(I,J,Z,3)
          BB9(J) = VBACK(I,J,Z,1)
          BB10(J)= VBACK(I,J,Z,2)
          BB11(J)= VBACK(I,J,Z,3)
        END DO
        AA1(I) = FILLIN(UP2,BB1,BB2,NCH)
        AA2(I) = FILLIN(UP2,BB1,BB3,NCH)
        AA3(I) = FILLIN(UP2,BB1,BB4,NCH)
        AA4(I) = FILLIN(UP2,BB1,BB5,NCH)
        AA5(I) = FILLIN(UP2,BB1,BB6,NCH)
        AA6(I) = FILLIN(UP2,BB1,BB7,NCH)
        AA7(I) = FILLIN(UP2,BB1,BB8,NCH)
        AA8(I) = FILLIN(UP2,BB1,BB9,NCH)
        AA9(I) = FILLIN(UP2,BB1,BB10,NCH)
        AA10(I)= FILLIN(UP2,BB1,BB11,NCH)
      END DO
      DEC_DU1P  = FILLIN(UP1,U11,AA1,NR)
      DEC_DU2P  = FILLIN(UP1,U11,AA2,NR)
      DET_DU1P  = FILLIN(UP1,U11,AA3,NR)
      DET_DU2P  = FILLIN(UP1,U11,AA4,NR)
     !Velocities in the cartesian system 
      VFACEP(1) = FILLIN(UP1,U11,AA5,NR)
      VFACEP(2) = FILLIN(UP1,U11,AA6,NR)
      VFACEP(3) = FILLIN(UP1,U11,AA7,NR)
      VBACKP(1) = FILLIN(UP1,U11,AA8,NR)
      VBACKP(2) = FILLIN(UP1,U11,AA9,NR)
      VBACKP(3) = FILLIN(UP1,U11,AA10,NR) 

     !Determinant of transformation Gij 
      G    =  G11P*(G22P*G33P-G23P*G32P)-G12P*(G21P*G33P-G23P*G31P)+G13P*(G21P*G32P-G22P*G31P)
      C11P = (G22P*G33P - G23P*G32P)/G 
      C21P = (G23P*G31P - G21P*G33P)/G
      C31P = (G21P*G32P - G22P*G31P)/G           
      C12P = (G32P*G13P - G33P*G12P)/G 
      C22P = (G33P*G11P - G31P*G13P)/G 
      C32P = (G31P*G12P - G32P*G11P)/G
      C13P = (G12P*G23P - G13P*G22P)/G 
      C23P = (G13P*G21P - G11P*G23P)/G 
      C33P = (G11P*G22P - G12P*G21P)/G 
     !Just to be sure we calculate the metric C at CP analyticaly
     !C    = C11P*(C22P*C33P-C23P*C32P)-C12P*(C21P*C33P-C23P*C31P)+C13P*(C21P*C32P-C22P*C31P) 
     !WRITE(*,*) 'C*G=',C*G

     !Undisturbed velocity expressed with the moving system's components (back)
      WBACKP(1)=VBACKP(1)*UNIT1(1)+VBACKP(2)*UNIT1(2)+VBACKP(3)*UNIT1(3) 
      WBACKP(2)=VBACKP(1)*UNIT2(1)+VBACKP(2)*UNIT2(2)+VBACKP(3)*UNIT2(3)
      WBACKP(3)=VBACKP(1)*UNIT3(1)+VBACKP(2)*UNIT3(2)+VBACKP(3)*UNIT3(3) 
     !Undisturbed velocity expressed with the moving system's components (face)
      WFACEP(1)=VFACEP(1)*UNIT1(1)+VFACEP(2)*UNIT1(2)+VFACEP(3)*UNIT1(3) 
      WFACEP(2)=VFACEP(1)*UNIT2(1)+VFACEP(2)*UNIT2(2)+VFACEP(3)*UNIT2(3)
      WFACEP(3)=VFACEP(1)*UNIT3(1)+VFACEP(2)*UNIT3(2)+VFACEP(3)*UNIT3(3) 

     !Velocities in the general curvilinear system 
     !Mean Velocity of face and back - Velocity on the reference surface for small thickness blades
     !The minus is due to the relative moving coordinate systems (forward velocity in one is backwards velocity in the other)
      Q(1) = -(C11P*(WBACKP(1)+WFACEP(1)) + C12P*(WBACKP(2)+WFACEP(2)) + C13P*(WBACKP(3)+WFACEP(3)))/2.D0 !>q1< 
      Q(2) = -(C21P*(WBACKP(1)+WFACEP(1)) + C22P*(WBACKP(2)+WFACEP(2)) + C23P*(WBACKP(3)+WFACEP(3)))/2.D0 !>q2<
      Q(3) = -(C31P*(WBACKP(1)+WFACEP(1)) + C32P*(WBACKP(2)+WFACEP(2)) + C33P*(WBACKP(3)+WFACEP(3)))/2.D0 !>q3<
     !Difference of Velocities on face and back
      Q(4) = -(C31P*(WBACKP(1)-WFACEP(1)) + C32P*(WBACKP(2)-WFACEP(2)) + C33P*(WBACKP(3)-WFACEP(3))) !<q3> 
      IF (N==N1+1) THEN    
     !Physical Component of Q^i on axis (3) at the L.E
      Q3_PHYSICAL(M,Z) = (GG13*Q(1) + GG33*Q(3))/SQRT(GG33)
     !Unit vector of \underline{x}^2 with components in the x^i system at the L.E
      E2(M,Z,1) = G12P/SQRT(GG22)
      E2(M,Z,2) = G22P/SQRT(GG22)
      E2(M,Z,3) = G32P/SQRT(GG22)
      END IF
     !Source and sinks on the blade (away from the L.E) 
      IF (N<N1+1) THEN
      S(M,N,Z) = (Q(1)*DET_DU1P + Q(2)*DET_DU2P - Q(4))*SQRT(GP/ALPHAP) 
      DPRE_S(M,N,Z) = Q(3)*SQRT(GP/ALPHAP)*S(M,N,Z) 
      END IF 
      DO J=1,4 ;QQ(M,N,Z,J)=Q(J) ;END DO

     !Normal vector of camber surface ,on the reference surface,written in cartesian coordinates
     !Notice that we dont divide by SQRT(ALPHA)
      NP(M,N,1) = (C11P*(-DEC_DU1P) + C21P*(-DEC_DU2P) + C31P*(1.D0))*SQRT(GP) 
      NP(M,N,2) = (C12P*(-DEC_DU1P) + C22P*(-DEC_DU2P) + C32P*(1.D0))*SQRT(GP) 
      NP(M,N,3) = (C13P*(-DEC_DU1P) + C23P*(-DEC_DU2P) + C33P*(1.D0))*SQRT(GP)                      
    
      END SUBROUTINE PROPERTIES_AT_CP


      
!     ******************************************************************
!     *                                                                *
!     *            SUBROUTINE LINEAR_SYSTEM_FORMATION                  *
!     *                                                                *
!     *       FOR THE CALCULATION OF ELEMENTS OF THE ARRAY             *
!     *                                                                *
!     ******************************************************************      
      
      SUBROUTINE LINEAR_SYSTEM_FORMATION

      INTEGER :: K,KZ,KZZ,II,JJ,L
      DOUBLE PRECISION :: RJJ

      IF (N<N1) THEN 
         CALL INDUCTION_FACTOR         
            
        !--------------------- NO - ENTRANCE BOUNDARY CONDITION ---------------------
         K = (M-1) + (N-1)*M0 + (Z-1)*M0*(P0+N0+1) 
         DO KZ=1,NBL
           KZZ = KZ + MOD(NBL-(Z-1),NBL)
           IF (KZZ>NBL) THEN 
           KZZ = KZZ - NBL
           END IF
           DO II=1,M0        
           L = II + (KZ-1)*M0*(N0+1)
           VT(K,L) = C31P*(TS_M0(M,N,KZZ,II,1) + 1*THT(M,N,KZZ,L,1)) &
                   + C32P*(TS_M0(M,N,KZZ,II,2) + 1*THT(M,N,KZZ,L,2)) &
                   + C33P*(TS_M0(M,N,KZZ,II,3) + 1*THT(M,N,KZZ,L,3)) !VT_SM0 
           END DO            
           DO II=1,M0 ;DO JJ=1,N0
           L = II + M0*JJ + (KZ-1)*M0*(N0+1)
           VT(K,L) = C31P*(TSS_MN(M,N,KZZ,II,JJ,1) + 1*THT(M,N,KZZ,L,1)) &
                   + C32P*(TSS_MN(M,N,KZZ,II,JJ,2) + 1*THT(M,N,KZZ,L,2)) &
                   + C33P*(TSS_MN(M,N,KZZ,II,JJ,3) + 1*THT(M,N,KZZ,L,3)) !VT_SSMN 
           END DO ;END DO
         END DO
       
        !Wake Induced Velocities
         IF ((IN==1).OR.(IN==2)) THEN
         WP(1)=0.D0 ;WP(2)=0.D0 ;WP(3)=0.D0 
         ELSE
         CALL WAKE_CONTRIBUTION
         END IF           
         B(K) = -(C31P*WP(1) + C32P*WP(2) + C33P*WP(3)) !to kanoniko thelei (-)
         B(K) = B(K) + (Q(1)*DEC_DU1P + Q(2)*DEC_DU2P - Q(3))
            
      ELSEIF (N==N1) THEN
   
        !-----------------TRAILING EDGE B.C - ZERO LOAD (BERNOULLI)-----------------------
         K  = (M-1) + (N-1)*M0 + (Z-1)*M0*(P0+N0+1)   
         II = M-1           
         IF (IN<=2) THEN  
         P = 1.D0/DT        !Zero laod at T.E (first order approximation of derivative)
         ELSE !(IN>=3)
         P = 3.D0/(2.D0*DT) !Zero laod at T.E (second order approximation of derivative)
         END IF            

         L = II + (Z-1)*M0*(N0+1)  !KZ=Z
         VT(K,L) = P*PI/2.D0 !VT_SM0
        !VT(K,L) = 0 !VT_SM0
         DO JJ=1,N0
         RJJ = DBLE(JJ)
         L = II + M0*JJ + (Z-1)*M0*(N0+1) !KZ=Z
         VT(K,L) = P*DBLE(1-(-1)**JJ)/RJJ + 2.D0*Q(2)*DBLE(JJ*(-1)**(JJ+1)) !VT_SSMN 
        !VT(K,L) = 2.D0*Q(2)*DBLE(JJ*(-1)**(JJ+1)) !VT_SSMN 
         END DO             

         IF (IN==1) THEN
            B(K) = 0.D0
         ELSE !(IN>=2)
            L = II + (Z-1)*M0*(N0+1) !KZ=Z
            IF (IN==2) THEN
            Q_SM0 =-A(L,IN-1)/DT
            ELSE !(IN>=3)
            Q_SM0 =(-4.D0*A(L,IN-1)+A(L,IN-2))/(2.D0*DT)
            END IF 
            B(K) = -(Q_SM0*PI/2.D0) 
           !B(K) = 0 
           
            DO JJ=1,N0
              RJJ = DBLE(JJ)
              L =  II + M0*JJ +(Z-1)*M0*(N0+1)                  
              IF (IN==2) THEN
              Q_SSMN =-A(L,IN-1)/DT
              ELSE
              Q_SSMN =(-4.D0*A(L,IN-1)+A(L,IN-2))/(2.D0*DT)
              END IF
              B(K) = B(K) - Q_SSMN*DBLE(1-(-1)**JJ)/RJJ 
             !B(K) = 0
            END DO
         END IF !IF (IN==1)
              
      ELSE !(N==0)-->(N==N1+1)
         WRITE(*,*) 'THERE ARE NO B.C AT THE L.E'
      END IF !IF (N<N1)      

      END SUBROUTINE LINEAR_SYSTEM_FORMATION     


      
!     ******************************************************************
!     *                                                                *
!     *              SUBROUTINE PRESSURE_DISTRIBUTION                  *
!     *                                                                *
!     *            FOR THE CALCULATION OF (P+-P-)*SQRT(a)              *
!     *                                                                *
!     ******************************************************************      

      SUBROUTINE PRESSURE_DISTRIBUTION

      INTEGER :: M,N,Z,II,JJ
      DOUBLE PRECISION :: RII,RJJ

     !---------- Calculation of (P+-P-)*a^(1/2) --------------
      DO Z=1,NBL           
      DO M=1,M1+1 
      DO N=1,N1+1 !Point N1+1 referes to the L.E
        DPRE(M,N,Z) = 0.D0
          
        IF ((M>1).AND.(M<M1+1)) THEN
        IF (IN==1) THEN
             DPRE(M,N,Z) = 0.D0
        ELSE !(IN>=2)
           YP1 = PI/DBLE(M1)*DBLE(M-1) 
          !Point N1+1 referes to the L.E
           IF (N==N1+1) THEN
             YP2 = 0.D0
           ELSE
             YP2 = PI/DBLE(N1)*DBLE(N)
           END IF
           DO J=1,4 ;Q(J)=QQ(M,N,Z,J) ;END DO

           IF (IN<=2) THEN 
           P = 1.D0/DT        !First order approximation of derivative
           ELSE !(IN>=3)
           P = 3.D0/(2.D0*DT) !Second order approximation of derivative
           END IF            

           DO II=1,M0
           RII = DBLE(II)  
           L = II + (Z-1)*M0*(N0+1) 
           IF (IN==2) THEN
           Q_SM0 =-A(L,IN-1)/DT
           ELSE !(IN>=3)
           Q_SM0 =(-4.D0*A(L,IN-1)+A(L,IN-2))/(2.D0*DT)
           END IF 
           DPRE(M,N,Z) = DPRE(M,N,Z) - (P*A(L,IN) + Q_SM0)*SIN(RII*YP1)*(YP2+SIN(YP2))/2.D0 &
                                 * 1.D0/4.D0 * SPAN * SIN(YP1)*SIN(YP2)     !DP_SM0 - Inertial pressure contribution        
           DO JJ=1,N0
              RJJ = DBLE(JJ)
              L =  II + M0*JJ +(Z-1)*M0*(N0+1)                  
              IF (IN==2) THEN
              Q_SSMN =-A(L,IN-1)/DT
              ELSE
              Q_SSMN =(-4.D0*A(L,IN-1)+A(L,IN-2))/(2.D0*DT)
              END IF
              DPRE(M,N,Z) = DPRE(M,N,Z) - (P*A(L,IN) + Q_SSMN)*SIN(RII*YP1)*(1-COS(RJJ*YP2))/RJJ &
                                    * 1.D0/4.D0 * SPAN * SIN(YP1)*SIN(YP2)  !DP_SSMN - Inertial pressure contribution
           END DO  
           END DO           

           DO II=1,M0
           RII = DBLE(II)  
           L = II + (Z-1)*M0*(N0+1) 
           DPRE(M,N,Z) = DPRE(M,N,Z) - Q(2)*A(L,IN)*SIN(RII*YP1)*(1.D0+COS(YP2)) &
                                 * 1.D0/4.D0 * SPAN * SIN(YP1)      !DP_SM0 - Dymamic pressure contribution         
           DO JJ=1,N0
              RJJ = DBLE(JJ)
              L =  II + M0*JJ +(Z-1)*M0*(N0+1)                  
              DPRE(M,N,Z) = DPRE(M,N,Z) - Q(2)*2.D0*A(L,IN)*SIN(RII*YP1)*SIN(RJJ*YP2) &
                                    * 1.D0/4.D0 * SPAN * SIN(YP1)   !DP_SSMN - Dymamic pressure contribution
           END DO        
           END DO      

           IF (N<N1+1) THEN
          !Contribution of sources-sinks to the pressure
           DPRE(M,N,Z) = DPRE(M,N,Z) - DPRE_S(M,N,Z)  
           END IF          
        END IF !IF (IN==1)
        END IF !IF ((M>1).AND.(M<M1+1)) THEN

        IF (MOD(M0,2)==0) THEN
        MID = (M0+2)/2
        ELSE
        MID = (M0+1)/2 +1
        END IF

        IF (IN==35) THEN
        IF ((M==MID).AND.(Z==1)) THEN
       !WRITE(5,*) YP2, DP(M,N,Z) 
        END IF
        END IF
      END DO !DO N=1,N1+1          
      END DO !DO M=1,M1+1 
      END DO !DO Z=1,NBL     


     !Calculation of leading edge suction force
      DO Z=1,NBL
      DO M=1,M1+1
        YP1 = PI/DBLE(M1)*DBLE(M-1) 
        GAMMA(M,Z) = 0.D0

       !For the Sinus Series, mu=0 at tip and hub            
        IF ((M>1).AND.(M<M1+1)) THEN             
        DO II=1,M0
          RII = DBLE(II)         
          L = II + (Z-1)*M0*(N0+1)
          GAMMA(M,Z) = GAMMA(M,Z) + A(L,IN)*SIN(RII*YP1)*PI/2.D0 !A_SM0
          DO JJ=1,N0
            RJJ = DBLE(JJ)
            L = II + M0*JJ + (Z-1)*M0*(N0+1)
            GAMMA(M,Z) = GAMMA(M,Z) + A(L,IN)*SIN(RII*YP1)*(1.D0-COS(RJJ*PI))/RJJ !A_SSMN                 
        END DO ;END DO
        
       !Leading Edge Suction Force
        LE_SF(M,Z) = Q3_PHYSICAL(M,Z)*GAMMA(M,Z) !-----------------------------------------------------------------******************************************Prosoxh sto prosimo
        LE_SF(M,Z) = LE_SF(M,Z)*0.5D0*SPAN*SIN(YP1)
        !WRITE(5,*)  Q3_PHYSICAL(M,Z),GAMMA(M,Z)*SPAN*0.5D0!(SIN(YP1))
       ELSE
        LE_SF(M,Z) = 0.D0
        END IF           
      END DO !M=1,M1+1
      END DO !Z=1,NBL
     
      END SUBROUTINE PRESSURE_DISTRIBUTION


       
      
!     ******************************************************************
!     *                                                                *
!     *                    SUBROUTINE LIFT_THRUST                      *
!     *                                                                *
!     *            FOR THE CALCULATION OF 2D-LIFT AND THRUST           *
!     *                                                                *
!     ******************************************************************

      SUBROUTINE LIFT_THRUST

      INTEGER :: M,N,Z,J
      
        DO J=1,3
        DO Z=1,NBL
        DO M=1,M1+1
        DO N=1,N1+1
          IF ((M==1).OR.(M==M1+1)) THEN
            !At tip and hub DP=0,and DP is not
            !defined in the program
            IF (N<N1+1) THEN
            AA11(N+1) = 0.D0
            ELSE !(N==N1+1)
            AA11(1)   = 0.D0
            END IF
          ELSE         
            IF (N<N1+1) THEN
            AA11(N+1) = DPRE(M,N,Z)* NP(M,N,J)
            ELSE !(N==N1+1)
            AA11(1)   = DPRE(M,N,Z)* NP(M,N,J)
            END IF
          END IF
        END DO !DO N=1,N1+1
        DL(M,Z) = FILON_METHOD(0.D0,PI,0,N1+1,AA11,0,1)
       
       !Correction of 2D-Lift by the leading edge suction force
        IF ((M>1).AND.(M<M1+1)) THEN
        DL(M,Z) = DL(M,Z) + LE_SF(M,Z)*E2(M,Z,J)
        END IF
        END DO !DO M=1,M1+1 
        END DO !DO Z=1,NBL    
 
        IF (MOD(M0,2)==0) THEN
        MID = (M0+2)/2
        ELSE
        MID = (M0+1)/2 +1
        END IF          
        IF (J==1)THEN
        WRITE(5,*) IN, DL(MID,1)/(0.5*SPAN)
        END IF

        DO Z=1,NBL
        DO M=1,M1+1
          AA12(M) = DL(M,Z)
        END DO !DO M=1,M1+1 
        LIFT(Z) = FILON_METHOD(0.D0,PI,0,M1+1,AA12,0,1)     
        END DO !DO Z=1,NBL    
            
        THRUST(IN,J) = 0.D0
        DO Z=1,NBL
        THRUST(IN,J) = THRUST(IN,J) + LIFT(Z)
        END DO
        
        IF (J==1) THEN
        !WRITE(5,*) IN, THRUST(IN,J) 
        END IF
        END DO !DO J=1,3

       !Thrust expressed with the inertial system's components   
        THRUST_I(IN,1) =THRUST(IN,1)*UNIT1(1) +THRUST(IN,2)*UNIT2(1) +THRUST(IN,3)*UNIT3(1)   
        THRUST_I(IN,2) =THRUST(IN,1)*UNIT1(2) +THRUST(IN,2)*UNIT2(2) +THRUST(IN,3)*UNIT3(2)
        THRUST_I(IN,3) =THRUST(IN,1)*UNIT1(3) +THRUST(IN,2)*UNIT2(3) +THRUST(IN,3)*UNIT3(3)  

        J=3
        !WRITE(5,*) IN, THRUST_I(IN,J)
      
      END SUBROUTINE LIFT_THRUST      
      
      

!     ******************************************************************
!     *                                                                *
!     *                 SUBROUTINE INDUCTION FACTOR                    *
!     *                                                                *
!     *  RETURNS THE INDUCTION FACTORS AND THT FOR THE SPECIFIC C.P    *
!     *                                                                *
!     ******************************************************************

      SUBROUTINE INDUCTION_FACTOR

      INTEGER KZ,II,JJ,J
      DOUBLE PRECISION :: RII,RJJ

      DO KZ=1,NBL    
        DELTA = 2.D0*PI*(DBLE(KZ)-1.D0)/DBLE(NBL) !The integration points, Q ,passes from all the blades            
        IF ((IN==1).AND.(Z==1)) THEN
       !Since the induction factors depend on the geometry only (and not time) we calculate them only the first time. !--------------------------------------------------------------------------------------Mporei na bgei eksw apo thn yporoutina
       !Since the propeller is symmetrical around its axis and the solution 1/r^3 of the laplacian is also symmetrical
       !the induction factors are the same for all the blades, so we calculate the only for the reference blade.  
            
       !Calculation of self induction factors BSJM,BCJM,CJ0,CJN      
        IF (KZ==1) THEN
       !Calculation of e^(jil)*Gi1*Gl2
        E(1) = G21P*G32P - G31P*G22P
        E(2) = G31P*G12P - G11P*G32P
        E(3) = G11P*G22P - G21P*G12P
        ELSE
        E(1)=0.D0 ;E(2)=0.D0 ;E(3)=0.D0 !self induction factor is zero for Q not on the reference blade
        END IF              

        DO J=1,3 
          DO II=1,M0
          RII = DBLE(II)
         !Calculation of BSJM_self
          EBS(II,J) = -8.D0*SIN(RII*YP1)*SQRT(A11P)/(ALPHAP*SPAN)*E(J)
          END DO 
         !Calculation of CJ0_self 
          EC0(J)= 4.D0*(YP2+SIN(YP2))*SQRT(A22P)/(ALPHAP*SPAN)*E(J) 
          DO JJ=1,N0
          RJJ = DBLE(J)
         !Calculation of CJN_self
          ECS(JJ,J) = 8.D0*(1.D0-COS(RJJ*YP2))*SQRT(A22P)/(ALPHAP*SPAN)*E(J)
          END DO 
        END DO

       !---------- Calculation of BSM0,BCM0,BSMN,BCMN,CSM0,CCM0,CSMN,CCMN -----------
        IF (KZ.NE.1) THEN
         !For KZ/=1 ,the self induction factors are zero so we ommit these terms
          CALL BIND (1,21,0.D0,PI,SS_S0,SS_SSN)
          DO J=1,3 ;DO II=1,M0 
            BSM0(II,J) = SS_S0(II,J)
            DO JJ=1,N0 
            BSSMN(II,JJ,J) = SS_SSN(II,JJ,J)
            END DO
          END DO ;END DO                
           CALL CIND (1,27,0.D0,PI,SS_S0,SS_SSN)
          DO J=1,3 ;DO II=1,M0
            CSM0(II,J) = SS_S0(II,J)
            DO JJ=1,N0 
            CSSMN(II,JJ,J) = SS_SSN(II,JJ,J)
            END DO
          END DO ;END DO                 
        END IF
              
       !Subroutine BIND   
       !Number of integration points for the outer integral, for BIND        
        F0  = PI/12.D0
        F1  = F0/2.D0
        IE0 = 13
        AB1 = YP2-F1
        IF (AB1<0.001D0) THEN
          F1  = YP2
          F0  = 2.D0*F1
          AB1 = 0.D0
        END IF
        AB2 = PI-YP2-F1
        IF (AB2<0.001D0) THEN
          F1  = PI - YP2
          F0  = 2.D0*F1
          AB2 = 0.D0
          AB1 = YP2 - F1 
        END IF
        IE1 = AINT(AB1/0.125D0)
        IE2 = AINT(AB2/0.125D0)
        IF (IE1.LT.501) THEN ; IE1=501 ;END IF
        IF (IE2.LT.501) THEN ; IE2=501 ;END IF
        IF (MOD(IE1,2).EQ.0) THEN ;IE1=IE1+1 ;END IF
        IF (MOD(IE2,2).EQ.2) THEN ;IE2=IE2+1 ;END IF

        CALL BIND (1,IE1,0.D0,AB1,SS_S01,SS_SSN1)          
        Y2A = YP2 - F1
        CALL BIND (0,IE0,Y2A ,F0 ,SS_S00,SS_SSN0)          
        Y2A = YP2 + F1
        CALL BIND (1,IE2,Y2A ,AB2,SS_S02,SS_SSN2)
        
        DO J=1,3 ;DO II=1,M0                
          BSM0(II,J) = SS_S01(II,J) + SS_S00(II,J) + SS_S02(II,J) + PI*EBS(II,J)/2.D0
          DO JJ=1,N0
          RJJ = DBLE(JJ)
          BSSMN(II,JJ,J) = SS_SSN1(II,JJ,J) + SS_SSN0(II,JJ,J) + SS_SSN2(II,JJ,J) - &
                           PI*EBS(II,J)*COS(RJJ*YP2)/SIN(YP2)
          END DO
        END DO ;END DO

       !Subroutine CSIND
       !Number of integration points for the outer integral, for CIND
        F0  = PI/12.D0
        F1  = F0/2.D0
        IE0 = 45
        AB1 = YP1 - F1
        IF (AB1<0.001D0) THEN
          F1  = YP1
          F0  = 2.D0*F1
          AB1 = 0.D0
        END IF
        AB2 = PI - YP1 - F1
        IF (AB2<0.001D0) THEN
          F1  = PI - YP1
          F0  = 2.D0*F1
          AB2 = 0.D0
          AB1 = YP1 - F1
        END IF
        IE1 = AINT(AB1/0.10D0)
        IE2 = AINT(AB2/0.10D0)
        IF (IE1<701) THEN ;IE1=701 ;END IF
        IF (IE2<701) THEN ;IE2=701 ;END IF
        IF (MOD(IE1,2)==0) THEN ;IE1=IE1+1 ;END IF
        IF (MOD(IE2,2)==0) THEN ;IE2=IE2+1 ;END IF

        CALL CIND (1,IE1,0.D0,AB1,SS_S01,SS_SSN1)
        Y1A = YP1 - F1
        CALL CIND (0,IE0,Y1A ,F0 ,SS_S00,SS_SSN0)
        Y1A = YP1 + F1
        CALL CIND (1,IE2,Y1A ,AB2,SS_S02,SS_SSN2)
        
        DO J=1,3 ;DO II=1,M0   
          RII = DBLE(II)                  
          CSM0(II,J) = SS_S01(II,J) + SS_S00(II,J) + SS_S02(II,J) + &
                       PI*EC0(J)*SIN(RII*YP1)/SIN(YP1)
          DO JJ=1,N0             
          CSSMN(II,JJ,J) = SS_SSN1(II,JJ,J) + SS_SSN0(II,JJ,J) + SS_SSN2(II,JJ,J) + &
                           PI*ECS(JJ,J)*SIN(RII*YP1)/SIN(YP1)
          END DO
        END DO ;END DO

       !Calculation of Mode Velocities 
        DO J=1,3 ;DO II=1,M0 
          RII = DBLE(II)
          TS_M0(M,N,KZ,II,J)=(0.5D0*SPAN*BSM0(II,J) -1*0.5D0*RII*CSM0(II,J))/(-4.D0*PI)   
          DO JJ=1,N0
          RJJ = DBLE(JJ)
          TSS_MN(M,N,KZ,II,JJ,J)=(0.5D0*SPAN*BSSMN(II,JJ,J) -1*0.5D0*RII/RJJ*CSSMN(II,JJ,J))/(-4.D0*PI)
          END DO
        END DO ;END DO             
        END IF !IF ((IN==1).AND.(Z==1)) THEN

        IF ((IN>=2).AND.(Z==1)) THEN
         !For IN==1 there is no wake or kutta strip
         !Contibution of the boundary line (Tip,Hub,TE-kutta strip) of KZ blade at control point P
          CALL TIP_HUB_TE(KZ)
          DO L=1,MNT 
          THT(M,N,KZ,L,:) = THTP(L,:) 
          END DO 
        ELSE
          DO L=1,MNT 
          THT(M,N,KZ,L,1) = 0.D0
          THT(M,N,KZ,L,2) = 0.D0
          THT(M,N,KZ,L,3) = 0.D0
          END DO
        END IF
      END DO !DO KZ=1,NBL      
       
      END SUBROUTINE INDUCTION_FACTOR      
      

!     ******************************************************************
!     *                                                                *
!     *                      SUBROUTINE BIND                           *
!     *                                                                *
!     *                FOR   THE CALCULATION OF BJM                    *
!     *                                                                *
!     ******************************************************************

      SUBROUTINE BIND (IX,IE,Y2A,AB,SS_S0,SS_SSN)
                     
      INTEGER ,INTENT(IN) :: IX,IE
      DOUBLE PRECISION ,INTENT(IN) :: Y2A,AB
      DOUBLE PRECISION ,INTENT(INOUT) :: SS_S0(M0,3),SS_SSN(M0,N0,3)    
     !Locals
      INTEGER :: L0,L,ND
      DOUBLE PRECISION :: EE,HC,HP,DY,DY1,DY2,DY0,DDY0,A0,B0,D0,R_R,R0,D
      DOUBLE PRECISION :: F(50201,3),FF(50201),FFF(IE),S_S(3),BS(3),R3,Y2B
      DOUBLE PRECISION ,DIMENSION(M0,IE,3) :: FBS0,FBSSN

      EE  = 0.0000001D0 !A very small quantity 
      HC  = AB/DBLE(IE-1) !Equal chordwise distance ,h, between integration points,Q

      DO L=1,IE
        Y2 = Y2A + HC*DBLE(L-1)
        U2 =-0.5D0*COS(Y2)
        DY = ABS(Y2-YP2)
        IF ((DY>0.005D0).OR.(DELTA>EE)) THEN
        !For (DELTA>0) C.P ,P,and integration point Q, belong on different blades
        !For (DY>0.005) C.P ,P,and integration point Q, are sufficiently far apart 
        !even on the same blade
          ! Number of inner intregral nodes
           IF (DELTA<=EE) THEN !DELTA=0
             DY1 = ABS(YP1-PI/2.D0)
             DY2 = ABS(YP2-PI/2.D0)
             DY0 = 1.2D0 - 2.D0*DY1/PI
             IF (DY<=DY0) THEN
               A0 = 40.D0 + 180.D0*DY1/PI - 500.D0*DY1**2/PI**2
               B0 = 0.9888D0 - 0.1876D0*DY1 + 0.1231D0*DY1**2
               D0 = 41.D0 + A0*(DY0/DY-1.D0)**B0
             ELSE
               D0 =41.D0
             END IF
             DDY0 = DY0 + 5.D0*DY2**2/PI**2
             IF (DY1<=DDY0) THEN
               R_R = 1.D0 + 0.39D0*DY2 - 1.025D0*DY2**2 + 1.249D0*DY2**3
             ELSE
               R_R = 1.D0
             END IF
             IF (DY1>=0.7685) THEN
               R0 = 0.585D0 + 0.54D0*DY1
             ELSE 
               R0 = 1.D0
             END IF
             D  = D0*R_R*R0
             ND = AINT(0.95D0*D)
             IF (MOD(ND,2)==0) THEN ;ND=ND+1 ;END IF
           ELSE
             ND = 33
           END IF ! IF (DELTA==0)
           IF (ND<401) THEN ;ND=401 ;END IF
           HP = PI/DBLE(ND-1) !Equal spanwise distance ,h, between integration points,Q

          !---- DEFINTION OF SPANWISE INTEGRANDS FOR CALCULATION OF BJM --------  
           DO I=1,ND
             Y1 = HP*DBLE(I-1)
             U1 = 0.5D0*(ROUT+RHUB) - 0.5D0*SPAN*COS(Y1)
            ! Calutlating geometric values at integration point,Q
             CHR = FILLIN(U1,U11,CHORD1,NR)  
             XS  = FILLIN(U1,U11,XS1 ,NR)
             SK  = FILLIN(U1,U11,SKEW1 ,NR)
             FI  = FILLIN(U1,U11,FI1 ,NR)
             DCHR= FILLIN(U1,U11,DCHORD1,NR)
             DXS = FILLIN(U1,U11,DXS1,NR)
             DSK = FILLIN(U1,U11,DSKEW1,NR)
             DFI = FILLIN(U1,U11,DFI1,NR)
              CALL GEOM(NBL,U1,U2,DELTA,CHR,XS,SK,FI,DCHR,DXS,DSK,DFI, &
             X1,X2,X3,ARG,G11,G21,G31,G12,G22,G32,G13,G23,G33)           
             R3= SQRT((X1-XP1)**2 + (X2-XP2)**2 + (X3-XP3)**2)**3
            !Calculation of sin(y1)*e^jil*Gi1*dx^l/r^3
             F(I,1) = SIN(Y1)*(G21*(X3-XP3) - G31*(X2-XP2))/R3
             F(I,2) = SIN(Y1)*(G31*(X1-XP1) - G11*(X3-XP3))/R3
             F(I,3) = SIN(Y1)*(G11*(X2-XP2) - G21*(X1-XP1))/R3
           END DO

           DO J=1,3 ;DO II=1,M0
             DO I=1,ND ;FF(I)=F(I,J) ;END DO
             S_S(J) = FILON_METHOD(DBLE(0),PI,II,ND,FF,1,2)
             IF (DELTA<=EE) THEN !DELTA==0
                BS(J) = S_S(J)*(COS(Y2)-COS(YP2))
                FBS0(II,L,J) = (BS(J)-EBS(II,J))/(COS(Y2)-COS(YP2))*(1+COS(Y2))/2.D0
                IF ((Y2>0.0001D0).AND.((PI-Y2)>0.0001D0)) THEN 
                ! Since for the no-entrance b.c the C.P isn't on the boundary ,if the 
                ! integration point Q isn't on the boundary then it could be Q=C.P  
                FBSSN(II,L,J)= (BS(J)/SIN(Y2)-EBS(II,J)/SIN(YP2))/(COS(Y2)-COS(YP2))*SIN(Y2)
                ELSE
                ! Since SIN(Y2)~0 ,then SELF(II,J)/SIN(YP2)*SIN(Y2)~0 ,so...
                FBSSN(II,L,J)= S_S(J)
                END IF
             ELSE
                FBS0(II,L,J) = S_S(J)*(1.D0+COS(Y2))/2.D0
                FBSSN(II,L,J)= S_S(J)
             END IF
           END DO ;END DO
        ELSE  
           L0 = L
        END IF  ! IF (DY>0.005D0).OR.(DELTA>EE)
      END DO !L=1,IE


      IF ((DELTA<EE).AND.(IX==0)) THEN  
      !This block could be under the L0=L but it doesn't matter because
      !for Y2<=0.005D0<F0=PI/12.D0 ,IX=0 (this works only for a single L0)    
        DO J=1,3 ;DO II=1,M0
          FBS0(II,L0,J) = 0.D0
          FBSSN(II,L0,J)= 0.D0
        END DO ;END DO
      END IF

      Y2B = Y2A+AB
      DO J=1,3 ;DO II=1,M0        
        DO L=1,IE ;FFF(L)= FBS0(II,L,J) ;END DO
        SS_S0(II,J) = FILON_METHOD(Y2A,Y2B,0,IE,FFF,IX,1)
        DO JJ=1,N0
           DO L=1,IE ;FFF(L)= FBSSN(II,L,J) ;END DO
           SS_SSN(II,JJ,J)= FILON_METHOD(Y2A,Y2B,JJ,IE,FFF,IX,2)
        END DO
      END DO ;END DO

      END SUBROUTINE BIND
      


!     ******************************************************************
!     *                                                                *
!     *                      SUBROUTINE CIND                           *
!     *                                                                *
!     *              FOR  THE CALCULATION OF  CJ0,CJN                  *
!     *                                                                *
!     ******************************************************************

      SUBROUTINE CIND (IX,IE,Y1A,AB,SS_S0,SS_SSN)

      INTEGER ,INTENT(IN) :: IX,IE
      DOUBLE PRECISION ,INTENT(IN) :: Y1A,AB
      DOUBLE PRECISION ,INTENT(INOUT) :: SS_S0(M0,3),SS_SSN(M0,N0,3) 
     !Locals
      INTEGER :: L0,L,ND
      DOUBLE PRECISION :: EE,HP,HC,DY,DY1,DY2,DY0,D00,D0,D ,R3,Y1B
      DOUBLE PRECISION ,DIMENSION(3) :: S_0,S_SN1,S_SN2,S_SN,C0,CSN 
      DOUBLE PRECISION :: F0(18700,3),FN(18700,3),FF(18700),FFF(IE)
      DOUBLE PRECISION :: FCS0(IE,3),FCSSN(N0,IE,3)     

      EE  = 0.0000001D0 !A very small quantity
      HP  = AB/DBLE(IE-1) !Equal spanwise distance ,h, between integration points,Q
 
      DO L=1,IE
        Y1 = Y1A + HP*DBLE(L-1)
        U1 = 0.5D0*(ROUT+RHUB) - 0.5D0*SPAN*COS(Y1)
               
       ! Calutlating geometric values at integration point,Q
        CHR = FILLIN(U1,U11,CHORD1,NR)  
        XS  = FILLIN(U1,U11,XS1 ,NR)
        SK  = FILLIN(U1,U11,SKEW1 ,NR)
        FI  = FILLIN(U1,U11,FI1 ,NR)
        DCHR= FILLIN(U1,U11,DCHORD1,NR)
        DXS = FILLIN(U1,U11,DXS1,NR)
        DSK = FILLIN(U1,U11,DSKEW1,NR)
        DFI = FILLIN(U1,U11,DFI1,NR)

        DY = ABS(Y1-YP1)
        IF ((DY>0.001D0).OR.(DELTA>EE)) THEN
       !For (DELTA>0) C.P ,P,and integration point Q, belong on different blades
       !For (DY>0.001) C.P ,P,and integration point Q, are sufficiently far apart 
       !even on the same blade
         ! Number of inner intregral nodes
          IF (DELTA<=EE) THEN !DELTA=0
            DY1 = ABS(YP1-PI/2.D0)
            DY2 = ABS(YP2-PI/2.D0)
            IF (DY<PI/20.D0) THEN
              D00 = 12.D0*(PI/20.D0/DY-1.D0)**1.12D0
              DY0 = PI/20.D0-0.076*DY2**1.678D0
              D0  = 21.D0+13.537D0*D00*(DY0**1.408D0)*EXP(0.105D0*DY2)
            ELSE 
              D0  = 21.D0
            END IF
            D  = D0*(1.D0+1.104D0*DY1-3.084D0*DY1**2 + 3.167D0*DY1**3)
            ND = AINT(D)
            IF (MOD(ND,2)==0) THEN ;ND=ND+1 ;END IF
          ELSE
            ND = 21
          END IF ! IF (DELTA=0)
          IF (ND<201) THEN ;ND=201 ;END IF
          HC = PI/DBLE(ND-1)  !Equal chordwise distance ,h, between integration points,Q
       
         ! DEFINITION OF CHORDWISE INTEGRANDS FOR CALCULATION OF CJ0,CJN
          DO J=1,ND
            Y2 = HC*DBLE(J-1)
            U2 =-0.5D0*COS(Y2)
             CALL GEOM(NBL,U1,U2,DELTA,CHR,XS,SK,FI,DCHR,DXS,DSK,DFI, &
            X1,X2,X3,ARG,G11,G21,G31,G12,G22,G32,G13,G23,G33)           
            R3= SQRT((X1-XP1)**2 + (X2-XP2)**2 + (X3-XP3)**2)**3
           !Calculation of sin(y2)*e^jil*Gi2*dx^l/r^3
            FN(J,1) = SIN(Y2)*(G22*(X3-XP3) - G32*(X2-XP2))/R3
            FN(J,2) = SIN(Y2)*(G32*(X1-XP1) - G12*(X3-XP3))/R3
            FN(J,3) = SIN(Y2)*(G12*(X2-XP2) - G22*(X1-XP1))/R3
            F0(J,1) = FN(J,1)*(Y2+SIN(Y2))/2.D0
            F0(J,2) = FN(J,2)*(Y2+SIN(Y2))/2.D0
            F0(J,3) = FN(J,3)*(Y2+SIN(Y2))/2.D0
          END DO


          DO J=1,3  
          DO I=1,ND ;FF(I) = F0(I,J) ;END DO
          S_0(J) = FILON_METHOD(DBLE(0),PI,0,ND,FF,1,1)
          DO I=1,ND ;FF(I) = FN(I,J) ;END DO
          S_SN1(J) = FILON_METHOD(DBLE(0),PI,0,ND,FF,1,1)   

          IF (DELTA<=EE) THEN ! DELTA=0
            C0(J) = S_0(J)*(COS(Y1)-COS(YP1))
            FCS0(L,J) = (C0(J)-EC0(J))/(COS(Y1)-COS(YP1))
          ELSE
            FCS0(L,J) = S_0(J)
          END IF  
           
          DO JJ=1,N0
            DO I=1,ND ;FF(I) = FN(I,J) ;END DO
            S_SN2(J) = FILON_METHOD(DBLE(0),PI,JJ,ND,FF,1,1)
            S_SN(J)  = S_SN1(J) - S_SN2(J)
            IF (DELTA<=EE) THEN !DELTA=0
              CSN(J) = S_SN(J)*(COS(Y1)-COS(YP1))
              FCSSN(JJ,L,J) = (CSN(J)-ECS(JJ,J))/(COS(Y1)-COS(YP1))
            ELSE
              FCSSN(JJ,L,J) = S_SN(J)
            END IF
          END DO
          END DO !DO J=1,3              
        ELSE
          L0=L
        END IF ! IF (DY>0.001).OR.(DELTA>EE)
      END DO !DO L=1,IE

      DO J=1,3
        IF ((DELTA<=EE).AND.(IX==0)) THEN
         !This block could be under the L0=L but it doesn't matter because
         !for Y1<=0.005<F0=PI/12. ,IX=0 (this works only for a single L0)  
          FCS0(L0,J) = 0.D0
          DO JJ=1,N0
          FCSSN(JJ,L0,J) = 0.D0        
          END DO 
        END IF 
        Y1B = Y1A + AB
        DO II=1,M0
          DO L=1,IE ;FFF(L)=FCS0(L,J) ;END DO
          SS_S0(II,J) = FILON_METHOD(Y1A,Y1B,II,IE,FFF,IX,1)
        END DO
        DO JJ=1,N0
          DO L=1,IE ;FFF(L)=FCSSN(JJ,L,J) ;END DO
          DO II=1,M0
            SS_SSN(II,JJ,J) = FILON_METHOD(Y1A,Y1B,II,IE,FFF,IX,1)
          END DO
        END DO
      END DO !DO J=1,3

      END SUBROUTINE CIND


      
!     ******************************************************************
!     *                                                                *
!     *                      FUNCTION    FILLIN                        *
!     *                                                                *
!     *                FOR   PARABOLIC   INTERPOLATION                 *
!     *                                                                *
!     ******************************************************************
    
	  FUNCTION FILLIN(X,AB,OR,N) RESULT(FILL)

      INTEGER :: N 
      DOUBLE PRECISION,INTENT(IN) ::  AB(N),OR(N)  
      DOUBLE PRECISION :: FILL ,X 
    
	  !Locals
      INTEGER :: I
      DOUBLE PRECISION :: X1,X2,X3 ,Y1,Y2,Y3 ,LN_X1,LN_X2,LN_X3
	  
      I=1;
      DO WHILE ((X>AB(I)).AND.(I<N))
          I = I+1 ;
      END DO

      IF (I==1) THEN !X<=AB(1)
          X1= AB(1) ;X2= AB(2) ;X3= AB(3);
          Y1= OR(1) ;Y2= OR(2) ;Y3= OR(3);

      ELSE IF (I==N) THEN !AB(N-1)<X<=AB(N) or X>AB(N)
          X1= AB(N-2) ;X2= AB(N-1) ;X3= AB(N);
          Y1= OR(N-2) ;Y2= OR(N-1) ;Y3= OR(N);
    
      ELSE !AB(I-1)<X<=AB(I) ,1<I<N,
          X1= AB(I-1) ;X2= AB(I) ;X3= AB(I+1);
          Y1= OR(I-1) ;Y2= OR(I) ;Y3= OR(I+1);
      END IF

      !2nd order lagrange polynomials
      LN_X1 = (X-X2)*(X-X3)/ ((X1-X2)*(X1-X3));
      LN_X2 = (X-X1)*(X-X3)/ ((X2-X1)*(X2-X3));
      LN_X3 = (X-X1)*(X-X2)/ ((X3-X1)*(X3-X2));

      !2nd order polynomial lagrange interpolation
      FILL = Y1*ln_X1 + Y2*ln_X2 + Y3*ln_X3 ;  
	  
      END FUNCTION FILLIN	 



!     ******************************************************************
!     *                                                                *
!     *                    SUBROUTINE GEOM                             *
!     *                                                                *
!     *   FOR THE CALCULATION OF METRIC OF BLADE REFERENCE SURFACE     *
!     *                                                                *
!     ******************************************************************

      SUBROUTINE GEOM(NBL,U1,U2,DLA,CHD,XS,SK,FI,DCHD,DXS,DSK,DFI,& 
                      X1,X2,X3,ARG,G11,G21,G31,G12,G22,G32,G13,G23,G33)

      INTEGER ,INTENT(IN) :: NBL
      DOUBLE PRECISION ,INTENT(IN):: U1,U2,DLA,CHD,XS,SK,FI,DCHD,DXS,DSK,DFI
      DOUBLE PRECISION ,INTENT(INOUT):: X1,X2,X3,ARG
      DOUBLE PRECISION ,INTENT(INOUT):: G11,G21,G31,G12,G22,G32,G13,G23,G33
      !Locals
      DOUBLE PRECISION :: GG

      !If we have only one blade we assume that we don't  
      !have a propeller but a wing
      IF (NBL==1) THEN
      !Cartesian Coordinates
      X1  = XS + CHD*SIN(FI)*U2
      X2  = U1
      X3  = U1*TAN(SK) + CHD*COS(FI)*U2
      !Components of Transformation matrix
      G11 = DXS + (DCHD*SIN(FI) + CHD*COS(FI)*DFI)*U2
      G21 = 1.D0
      G31 = TAN(SK) + (U1/(COS(SK))**2)*DSK + (DCHD*COS(FI) - CHD*SIN(FI)*DFI)*U2
      G12 = CHD*SIN(FI)
      G22 = 0.D0
      G32 = CHD*COS(FI)
      G13 =-CHD*COS(FI)
      G23 = 0.D0
      G33 = CHD*SIN(FI) 

      ELSE  
      !Cartesian Coordinates
      ARG = SK + CHD*COS(FI)*U2/U1 + DLA
      X1  = XS + CHD*SIN(FI)*U2
      X2  = U1*COS(ARG)
      X3  = U1*SIN(ARG)
      !Components of Transformation matrix
      G11 = DXS + (DCHD*SIN(FI) + CHD*COS(FI)*DFI)*U2
      GG  = DSK + (DCHD*COS(FI) - CHD*SIN(FI)*DFI - CHD*COS(FI)/U1)*U2/U1
      G21 = COS(ARG) - U1*SIN(ARG)*GG
      G31 = SIN(ARG) + U1*COS(ARG)*GG
      G12 = CHD*SIN(FI)
      G22 =-CHD*COS(FI)*SIN(ARG)
      G32 = CHD*COS(FI)*COS(ARG)
      G13 =-CHD*COS(FI)
      G23 =-CHD*SIN(FI)*SIN(ARG)
      G33 = CHD*SIN(FI)*COS(ARG)
      END IF
     
      END SUBROUTINE GEOM 


!     ******************************************************************
!     *                                                                *
!     *                   FUNCTION  FILON  METHOD                      *
!     *                                                                *
!     *   FOR INTEGRATING FUNCTIONS CONTAINING TRIGONOMETRIC TERMS     *
!     *                                                                *
!     ******************************************************************      

      FUNCTION FILON_METHOD(a,b,p,N,f,filon,COS_SIN) RESULT(FILON_METH)

      INTEGER :: N ,p ,filon,COS_SIN
      DOUBLE PRECISION :: FILON_METH , a,b,f(N)   !<------------------------------------------ Prosexe tis diastaseis na mhn einai metablhtes
	  
	  !Locals
      INTEGER :: I
      DOUBLE PRECISION :: C_K,S_K  ,C_2K_1,C_2K  ,S_2K_1,S_2K
      DOUBLE PRECISION :: alpha,beta,gamma ,h,t ,x,EE
      DOUBLE PRECISION :: I1,I2,J1,J2
	  
      EE = 0.0001D0 !small quantity
      	  
      h = (b-a)/DBLE(N-1)
      t = p*h  !t= theta
	  
      IF (filon==0) THEN
         ! ---------------- 1st order Filon ---------------
         IF (t>EE) THEN
             alpha = t**(-1) - t**(-2.D0)*sin(t) 
             beta  = 2.D0*t**(-2)*(1.D0-cos(t)) 
         ELSE
             alpha = t/6.D0 - t**3/120.D0 + t**5/5040.D0
             beta  = 1.D0 - t**2/12.D0 + t**4/360.D0 - t**6/20160.D0 
         END IF		 
         ! ----- I = integral(a,b)(f(x)cos(px))dx ---------
         IF (COS_SIN==1) THEN
             C_K = 0
             DO I=2,N-1
                 x = a + h*DBLE(I-1)
                 C_K = C_K + f(I)*cos(p*x) 
             END DO
             C_K = C_K + 0.5D0*(f(1)*cos(p*a) + f(N)*cos(p*b))
             
             I1 = h*(alpha*(f(N)*sin(p*b) - f(1)*sin(p*a)) + beta*C_K)
             FILON_METH  = I1
            !write(*,*) 'I1=',I1
         END IF
         ! ----- J = integral(a,b)(f(x)sin(px))dx ---------
         IF (COS_SIN==2) THEN
             S_K = 0
             DO I=2,N-1
                 x = a + h*DBLE(I-1)
                 S_K = S_K + f(I)*sin(p*x) 
             END DO
             S_K = S_K + 0.5D0*(f(1)*sin(p*a) + f(N)*sin(p*b))

             J1 = h*(-alpha*(f(N)*cos(p*b) - f(1)*cos(p*a)) + beta*S_K)
             FILON_METH  = J1
            !write(*,*) 'J1=',J1
         END IF
      END IF

      IF (filon==1) THEN
         ! ---------------- 2nd order Filon ---------------
         IF (t>EE) THEN
             alpha = (t**2 + t*sin(t)*cos(t) - 2.D0*(sin(t))**2 ) /t**3 
             beta  = 2.D0*( t*(1+(cos(t))**2) - 2.D0*sin(t)*cos(t) )/t**3 
             gamma = 4.D0*(sin(t) - t*cos(t)) / t**3 
         ELSE
             alpha = 2.D0*t**3/45.D0 - 2.D0*t**5/315.D0 + 2.D0*t**7/4725.D0 - 8.D0*t**9/467775.D0 
             beta  = 2.D0/3.D0 + 2.D0*t**2/15.D0 - 4.D0*t**4/105.D0 + 2.D0*t**6/567.D0 
             gamma = 4.D0/3.D0 - 2.D0*t**2/15.D0 + t**4/210.D0 - t**6/11340.D0 
         END IF
         ! ----- I = integral(a,b)(f(x)cos(px))dx ---------
         IF (COS_SIN==1) THEN
             C_2K_1 = 0
             DO I=3,N-2,2
                 x = a + h*DBLE(I-1)
                 C_2K_1 = C_2K_1 + f(I)*cos(p*x) 
             END DO
             C_2K_1 = C_2K_1 + 0.5D0*(f(1)*cos(p*a) + f(N)*cos(p*b)) 

             C_2K = 0
             DO I=2,N-1,2
                 x = a + h*DBLE(I-1)
                 C_2K = C_2K + f(I)*cos(p*x)
             END DO
             I2 = h*(a*(f(N)*sin(p*b) - f(1)*sin(p*a)) + beta*C_2K_1 + gamma*C_2K)
             FILON_METH  = I2
             !write(*,*) 'I2=',I2
         END IF		 
         ! ----- J = integral(a,b)(f(x)sin(px))dx --------- 
         IF (COS_SIN==2) THEN
             S_2K_1 = 0
             DO I=3,N-2,2
                 x = a + h*DBLE(I-1)
                 S_2K_1 = S_2K_1 + f(I)*sin(p*x) 
             END DO
             S_2K_1 = S_2K_1 + 0.5D0*(f(1)*sin(p*a) + f(N)*sin(p*b)) 

             S_2K = 0
             DO I=2,N-1,2
                 x = a + h*DBLE(I-1)
                 S_2K = S_2K + f(I)*sin(p*x);
             END DO
             J2 = h*(-a*(f(N)*cos(p*b) - f(1)*cos(p*a)) + beta*S_2K_1 + gamma*S_2K)
             FILON_METH  = J2
             !write(*,*) 'J2=',J2
         END IF
      END IF
      	  
      END FUNCTION FILON_METHOD


      
!     ******************************************************************
!     *                                                                *
!     *                   SUBROUTINE TIP_HUB_TE                        *
!     *                                                                *
!     *      FOR THE CALCULATION OF TIP_HUB_TE INDUCED VELOCITIES      *
!     *                                                                *
!     ******************************************************************   

      SUBROUTINE TIP_HUB_TE(KZ)

      INTEGER,INTENT(IN) :: KZ
      !Locals
      INTEGER :: M,N,II,JJ,L,I,M1,N1
      DOUBLE PRECISION :: SMALL,FACTOR ,R,CC
      DOUBLE PRECISION :: H_AB,H_BC,H_CD,H_DA  ,A_AB,B_AB,B_BC,C_BC,C_CD,D_CD,D_DA,A_DA
      DOUBLE PRECISION ,DIMENSION(3) :: A,B,C,D,P  ,AB,BC,CD,DA  ,PA,PB,PC,PD
      DOUBLE PRECISION ,DIMENSION(3) :: PAxPB,PBxPC,PCxPD,PDxPA  ,E_AB,E_BC,E_CD,E_DA
      DOUBLE PRECISION ,DIMENSION(CPW,P0+N0+1,3) :: Q ,EP
      DOUBLE PRECISION ,DIMENSION(CPW,P0+N0+1,MNT_MAX) :: MU

      SMALL = 0.001D0
      P(1)=XP1 ;P(2)=XP2 ;P(3)=XP3
      I  = IN-1
      M1 = M0+1        !M1+1 C.P in the y1 direction including those on tip and hub   
      N1 = P0+N0+1     !N1+1 C.P in the y2 direction including those on L.E and T.E                 
     !Viscous Filter parameters
      R = 0.10D0 ![m] Viscous range
      CC = 0.69314718

      DO M=1,CPW
         !For M=1-->YP1=0 (hub) and M=CPW-->YP1=PI (tip)
          Y1 = PI/DBLE(CPW-1)*DBLE(M-1) 
          U1 = 0.5D0*(ROUT+RHUB) - 0.5D0*(ROUT-RHUB)*COS(Y1)       
         !Calutlating geometric values at integration point,Q
          CHR = FILLIN(U1,U11,CHORD1,NR)  
          XS  = FILLIN(U1,U11,XS1 ,NR)
          SK  = FILLIN(U1,U11,SKEW1 ,NR)
          FI  = FILLIN(U1,U11,FI1 ,NR)
          DCHR= FILLIN(U1,U11,DCHORD1,NR)
          DXS = FILLIN(U1,U11,DXS1,NR)
          DSK = FILLIN(U1,U11,DSKEW1,NR)
          DFI = FILLIN(U1,U11,DFI1,NR)      
      DO N=1,N1
         IF ((M==1).OR.(M==CPW).OR.(N==N1)) THEN
          !For N=0-->YP2=0 (L.E) and N=N1-->YP1=PI (T.E)
          Y2  = PI/DBLE(N1)*DBLE(N)     
          U2  =-0.5D0*COS(Y2) 
          DELTA = 2.D0*PI*(DBLE(KZ)-1.D0)/DBLE(NBL)  !The C.P belong on the reference blade    
           CALL GEOM(NBL,U1,U2,DELTA,CHR,XS,SK,FI,DCHR,DXS,DSK,DFI, &
          X1,X2,X3,ARG,G11,G21,G31,G12,G22,G32,G13,G23,G33)
          Q(M,N,1)=X1 ;Q(M,N,2)=X2 ;Q(M,N,3)=X3

          DO II=1,M0
             RII = DBLE(II)
             L = II + (KZ-1)*M0*(N0+1)
             MU(M,N,L) = SIN(RII*Y1)*(Y2+SIN(Y2))/2.D0 !MU_SM0 
             THTP(L,:) = 0.D0 !THTP_SM0
             DO JJ=1,N0
               RJJ = DBLE(JJ)
               L = II + M0*JJ + (KZ-1)*M0*(N0+1)
               MU(M,N,L) = SIN(RII*Y1)*(1.D0-COS(RJJ*Y2))/RJJ !MU_SSMN  
               THTP(L,:) = 0.D0 !THTP_SSMN  
          END DO ;END DO 

          END IF !IF ((M==1).OR.(M==CPW).OR.(N==N1)) THEN      
      END DO !DO N=1,N1
      END DO !DO M=1,CPW

 
      DO M=1,CPW
      DO N=1,N1
        IF ((M==CPW).AND.(N<N1)) THEN
        A(1)= Q(M,N+1,1) ;A(2)= Q(M,N+1,2) ;A(3)= Q(M,N+1,3)     !  L.E---B--<--A---T.E
        B(1)= Q(M,N,1)   ;B(2)= Q(M,N,2)   ;B(3)= Q(M,N,3)
        END IF   
        IF ((M==1).AND.(N<N1)) THEN
        A(1)= Q(M,N,1)   ;A(2)= Q(M,N,2)   ;A(3)= Q(M,N,3)       !  L.E---A-->--B---T.E
        B(1)= Q(M,N+1,1) ;B(2)= Q(M,N+1,2) ;B(3)= Q(M,N+1,3)      
        END IF
        
        IF (((M==CPW).OR.(M==1)).AND.(N<N1)) THEN
       !Vectors of line segments on the boundary
        DO J=1,3
         AB(J) = B(J) - A(J)
        END DO
       !Unit vectors of line segments on the boundary 
        DO J=1,3
         AB(J) = (B(J) - A(J))/NORM(AB)
        END DO    
        H_AB = POINT_TO_LINE_DISTANCE(P,A,AB)
        DO J=1,3
         PA(J) = A(J) - P(J)  
         PB(J) = B(J) - P(J)
        END DO
        PAxPB(:) = CROSS(PA,PB)   
       !Angles
        A_AB = VECTOR_ANGLE( AB,-PA)   
        B_AB = VECTOR_ANGLE(-AB,-PB)
        DO J=1,3
         IF (ABS(COS(A_AB)+COS(B_AB))<SMALL) THEN
         E_AB(J) = 0.D0
         ELSE
         E_AB(J) =-(COS(A_AB)+COS(B_AB))*PAxPB(J)/NORM(PAxPB)/(4*PI*H_AB)
         END IF
        END DO     
        DO J=1,3
         EP(M,N,J) = E_AB(J) 
        END DO
        END IF !IF (((M==CPW).OR.(M==1)).AND.(N<N1)) THEN
        

        IF ((M<CPW).AND.(N==N1)) THEN                                                    
        A(1)= Q(M+1,N,1)        ;A(2)= Q(M+1,N,2)        ;A(3)= Q(M+1,N,3)                
        B(1)= Q(M,N,1)          ;B(2)= Q(M,N,2)          ;B(3)= Q(M,N,3)                   !  A--<--D
        IF (IN>=2) THEN                                                                    !  |     |
        C(1)= WAKEP(M,KZ,I,1)   ;C(2)= WAKEP(M,KZ,I,2)   ;C(3)= WAKEP(M,KZ,I,3)            !  |     |
        D(1)= WAKEP(M+1,KZ,I,1) ;D(2)= WAKEP(M+1,KZ,I,2) ;D(3)= WAKEP(M+1,KZ,I,3)          !  B-->--C
        END IF
       !Vectors of line segments on the boundary
        DO J=1,3
         AB(J) = B(J) - A(J)
         IF (IN>=2) THEN
         BC(J) = C(J) - B(J)
         CD(J) = D(J) - C(J)
         DA(J) = A(J) - D(J)
         END IF
        END DO
       !Unit vectors of line segments on the boundary 
        DO J=1,3
         AB(J) = (B(J) - A(J))/NORM(AB)
         IF (IN>=2) THEN
         BC(J) = (C(J) - B(J))/NORM(BC) 
         CD(J) = (D(J) - C(J))/NORM(CD)
         DA(J) = (A(J) - D(J))/NORM(DA)
         END IF
        END DO    
        H_AB = POINT_TO_LINE_DISTANCE(P,A,AB)
        IF (IN>=2) THEN
        H_BC = POINT_TO_LINE_DISTANCE(P,B,BC)
        H_CD = POINT_TO_LINE_DISTANCE(P,C,CD)
        H_DA = POINT_TO_LINE_DISTANCE(P,D,DA)
        END IF
        DO J=1,3
         PA(J) = A(J) - P(J)
         PB(J) = B(J) - P(J)
         IF (IN>=2) THEN
         PC(J) = C(J) - P(J)
         PD(J) = D(J) - P(J)
         END IF
        END DO
        PAxPB(:) = CROSS(PA,PB)
        IF (IN>=2) THEN
        PBxPC(:) = CROSS(PB,PC)
        PCxPD(:) = CROSS(PC,PD)
        PDxPA(:) = CROSS(PD,PA)
        END IF
       !Angles [RAD]
        A_AB = VECTOR_ANGLE( AB,-PA)   
        B_AB = VECTOR_ANGLE(-AB,-PB)
        IF (IN>=2) THEN
        B_BC = VECTOR_ANGLE( BC,-PB)
        C_BC = VECTOR_ANGLE(-BC,-PC)
        C_CD = VECTOR_ANGLE( CD,-PC)
        D_CD = VECTOR_ANGLE(-CD,-PD)
        D_DA = VECTOR_ANGLE( DA,-PD)
        A_DA = VECTOR_ANGLE(-DA,-PA)
        END IF
        DO J=1,3
          IF (ABS(COS(A_AB)+COS(B_AB))<SMALL) THEN
          E_AB(J) = 0.D0  
          ELSE
          E_AB(J) =-(COS(A_AB)+COS(B_AB))*PAxPB(J)/NORM(PAxPB)/(4*PI*H_AB)!*(1-EXP(-CC*(H_AB/R)**2)) 
          END IF
          IF (IN>=2) THEN
          IF (ABS(COS(B_BC)+COS(C_BC))<SMALL) THEN
          E_BC(J) = 0.D0  
          ELSE
          E_BC(J) =-(COS(B_BC)+COS(C_BC))*PBxPC(J)/NORM(PBxPC)/(4*PI*H_BC)!*(1-EXP(-CC*(H_BC/R)**2)) 
          END IF
          IF (ABS(COS(C_CD)+COS(D_CD))<SMALL) THEN
          E_CD(J) = 0.D0
          ELSE
          E_CD(J) =-(COS(C_CD)+COS(D_CD))*PCxPD(J)/NORM(PCxPD)/(4*PI*H_CD)!*(1-EXP(-CC*(H_CD/R)**2)) 
          END IF
          IF (ABS(COS(D_DA)+COS(A_DA))<SMALL) THEN
          E_DA(J) = 0.D0
          ELSE  
          E_DA(J) =-(COS(D_DA)+COS(A_DA))*PDxPA(J)/NORM(PDxPA)/(4*PI*H_DA)!*(1-EXP(-CC*(H_DA/R)**2)) 
          END IF
          END IF
        END DO
        DO J=1,3
         IF (IN==1) THEN
         EP(M,N,J) = E_AB(J)   
         ELSE !(IN>=2)
         EP(M,N,J) =-E_AB(J) + E_BC(J) + E_CD(J) + E_DA(J)
         END IF
        END DO        
        END IF !IF ((M<CPW).AND.(N==N1))
                   
      END DO !DO N=1,N1
      END DO !DO M=1,CPW  


      DO J=1,3
      DO M=1,CPW
      DO N=1,N1
        
      IF ((M==1).AND.(N<N1)) THEN
        DO II=1,M0
        L = II + (KZ-1)*M0*(N0+1)
        THTP(L,J)= THTP(L,J) + 0*EP(M,N,J)*(MU(M,N+1,L)+MU(M,N,L))/2.D0 !THTP_SM0               
        DO JJ=1,N0
          L = II + M0*JJ + (KZ-1)*M0*(N0+1)
          THTP(L,J)= THTP(L,J) + 0*EP(M,N,J)*(MU(M,N+1,L)+MU(M,N,L))/2.D0 !THTP_SSMN 
        END DO ;END DO   
      END IF  
      IF ((M==CPW).AND.(N<N1)) THEN
        DO II=1,M0
          L = II + (KZ-1)*M0*(N0+1)
          THTP(L,J)= THTP(L,J) + 0*EP(M,N,J)*(MU(M,N+1,L)+MU(M,N,L))/2.D0 !THTP_SM0
          DO JJ=1,N0
            L = II + M0*JJ + (KZ-1)*M0*(N0+1)
            THTP(L,J)= THTP(L,J) + 0*EP(M,N,J)*(MU(M,N+1,L)+MU(M,N,L))/2.D0 !THTP_SSMN
        END DO ;END DO             
      END IF
      IF ((M<CPW).AND.(N==N1)) THEN
        IF (IN==1) THEN
          FACTOR = 2.D0
        ELSE !(IN>=2) 
          FACTOR = 4.D0
        END IF
        DO II=1,M0
          L = II + (KZ-1)*M0*(N0+1)
          THTP(L,J)= THTP(L,J) + EP(M,N,J)*(MU(M+1,N,L)+MU(M,N,L))/FACTOR !THTP_SM0
          DO JJ=1,N0
            L = II + M0*JJ + (KZ-1)*M0*(N0+1)
            THTP(L,J)= THTP(L,J) + EP(M,N,J)*(MU(M+1,N,L)+MU(M,N,L))/FACTOR !THTP_SSMN
        END DO ;END DO             
      END IF

      END DO !DO N=1,N1
      END DO !DO M=1,CPW 
      END DO !DO J=1,3 
      
      END SUBROUTINE TIP_HUB_TE

      

!     ******************************************************************
!     *                                                                *
!     *               SUBROUTINE WAKE_CONTRIBUTION                     *
!     *                                                                *
!     *        FOR THE CALCULATION OF WAKE INDUCED VELOCITIES          *
!     *                                                                *
!     ******************************************************************   

      SUBROUTINE WAKE_CONTRIBUTION                                                  

      !Locals
      INTEGER :: I,J,Z,M
      DOUBLE PRECISION :: MU ,H_AB,H_BC,H_CD,H_DA ,SMALL ,R,CC
      DOUBLE PRECISION :: A_AB,B_AB,B_BC,C_BC,C_CD,D_CD,D_DA,A_DA
      DOUBLE PRECISION ,DIMENSION(3) :: A,B,C,D,P  ,AB,BC,CD,DA ,PA,PB,PC,PD
      DOUBLE PRECISION ,DIMENSION(3) :: PAxPB,PBxPC,PCxPD,PDxPA ,W_AB,W_BC,W_CD,W_DA

      SMALL = 0.001D0
      P(1)=XP1 ;P(2)=XP2 ;P(3)=XP3
      WP(1)=0.D0 ;WP(2)=0.D0 ;WP(3)=0.D0
     !Viscous Filter parameters
      R = 0.10D0 ![m] Viscous range
      CC = 0.69314718
   
      !R = R - DBLE(IN)/60

      DO I=1,IN-1 !IN>=3 
      DO Z=1,NBL    
      DO M=1,CPW-1 
      
      !Wake Boundary element with constant dipole  
      !intesity equal to the mean value of its nodes  
      IF (I<IN-1) THEN
      MU = (M_WAKE(M+1,Z,I+1)+M_WAKE(M,Z,I+1)+M_WAKE(M,Z,I)+M_WAKE(M+1,Z,I))/4.D0
      ELSE !I=IN-1
      MU = (M_WAKE(M,Z,I)+M_WAKE(M+1,Z,I))/4.D0
      END IF
   
      !A,B,C,D are the position vectors for the nodes of the wake element
      A(1)= WAKEP(M+1,Z,I+1,1) ;A(2)= WAKEP(M+1,Z,I+1,2) ;A(3)= WAKEP(M+1,Z,I+1,3)    !   A--<--D
      B(1)= WAKEP(M,Z,I+1,1)   ;B(2)= WAKEP(M,Z,I+1,2)   ;B(3)= WAKEP(M,Z,I+1,3)      !   |     |
      C(1)= WAKEP(M,Z,I,1)     ;C(2)= WAKEP(M,Z,I,2)     ;C(3)= WAKEP(M,Z,I,3)        !   |     |
      D(1)= WAKEP(M+1,Z,I,1)   ;D(2)= WAKEP(M+1,Z,I,2)   ;D(3)= WAKEP(M+1,Z,I,3)      !   B-->--C

      !Vectors of line segments on the boundary of the wake element
      DO J=1,3
      AB(J) = B(J) - A(J)
      BC(J) = C(J) - B(J)
      CD(J) = D(J) - C(J)
      DA(J) = A(J) - D(J)
      END DO
      !Unit vectors of line segments on the boundary of the wake element
      DO J=1,3
      AB(J) = (B(J) - A(J))/NORM(AB)
      BC(J) = (C(J) - B(J))/NORM(BC) 
      CD(J) = (D(J) - C(J))/NORM(CD)
      DA(J) = (A(J) - D(J))/NORM(DA)
      END DO    
      H_AB = POINT_TO_LINE_DISTANCE(P,A,AB)
      H_BC = POINT_TO_LINE_DISTANCE(P,B,BC)
      H_CD = POINT_TO_LINE_DISTANCE(P,C,CD)
      H_DA = POINT_TO_LINE_DISTANCE(P,D,DA)

      DO J=1,3
      PA(J) = A(J) - P(J)  
      PB(J) = B(J) - P(J)
      PC(J) = C(J) - P(J)
      PD(J) = D(J) - P(J)
      END DO
      PAxPB(:) = CROSS(PA,PB) 
      PBxPC(:) = CROSS(PB,PC)
      PCxPD(:) = CROSS(PC,PD)
      PDxPA(:) = CROSS(PD,PA)

      !Angles [RAD]
      A_AB = VECTOR_ANGLE( AB,-PA)  
      B_AB = VECTOR_ANGLE(-AB,-PB)
      B_BC = VECTOR_ANGLE( BC,-PB)
      C_BC = VECTOR_ANGLE(-BC,-PC)
      C_CD = VECTOR_ANGLE( CD,-PC)
      D_CD = VECTOR_ANGLE(-CD,-PD)
      D_DA = VECTOR_ANGLE( DA,-PD)
      A_DA = VECTOR_ANGLE(-DA,-PA)
 
      DO J=1,3
       IF (ABS(COS(A_AB)+COS(B_AB))<SMALL) THEN
       W_AB(J) = 0.D0
       ELSE   
       W_AB(J) =-(COS(A_AB)+COS(B_AB))*PAxPB(J)/NORM(PAxPB)/(4*PI*H_AB)!*(1-EXP(-CC*(H_AB/R)**2))  
       END IF
       IF (ABS(COS(B_BC)+COS(C_BC))<SMALL) THEN
       W_BC(J) = 0.D0  
       ELSE  
       W_BC(J) =-(COS(B_BC)+COS(C_BC))*PBxPC(J)/NORM(PBxPC)/(4*PI*H_BC)!*(1-EXP(-CC*(H_BC/R)**2)) 
       END IF
       IF (ABS(COS(C_CD)+COS(D_CD))<SMALL) THEN
       W_CD(J) = 0.D0
       ELSE 
       W_CD(J) =-(COS(C_CD)+COS(D_CD))*PCxPD(J)/NORM(PCxPD)/(4*PI*H_CD)!*(1-EXP(-CC*(H_CD/R)**2)) 
       END IF
       IF (ABS(COS(D_DA)+COS(A_DA))<SMALL) THEN
       W_DA(J) = 0.D0  
       ELSE 
       W_DA(J) =-(COS(D_DA)+COS(A_DA))*PDxPA(J)/NORM(PDxPA)/(4*PI*H_DA)!*(1-EXP(-CC*(H_DA/R)**2)) 
       END IF
      END DO  
      DO J=1,3
      WP(J)   = WP(J) + MU*(W_AB(J) + W_BC(J) + W_CD(J) + W_DA(J))
      END DO
      
      END DO !DO M=1,CPW-1
      END DO !DO Z=1,NBL
      END DO !DO I=1,IN-1

      END SUBROUTINE WAKE_CONTRIBUTION

      
!----------------------------------------------------------------------      
	  FUNCTION NORM(V) RESULT(E_NORM)
      DOUBLE PRECISION ,INTENT(IN) :: V(3)
      DOUBLE PRECISION :: E_NORM
      E_NORM = SQRT(V(1)**2 + V(2)**2 +V(3)**2)
      END FUNCTION NORM

!----------------------------------------------------------------------
	  FUNCTION CROSS(V1,V2) RESULT(E_CROSS)
      DOUBLE PRECISION ,INTENT(IN):: V1(3),V2(3)
      DOUBLE PRECISION :: E_CROSS(3)
      E_CROSS(1) = V1(2)*V2(3) - V1(3)*V2(2)
      E_CROSS(2) = V1(3)*V2(1) - V1(1)*V2(3)
      E_CROSS(3) = V1(1)*V2(2) - V1(2)*V2(1)
      END FUNCTION CROSS

!---------------------------------------------------------------------
      FUNCTION POINT_TO_LINE_DISTANCE(P,A,E1) RESULT(R)
      DOUBLE PRECISION ,INTENT(IN) :: P(3) !Position vector at point P
      DOUBLE PRECISION ,INTENT(IN) :: A(3),E1(3) !define axis position by <point> and <direction unit vector>
      DOUBLE PRECISION :: R
      !Locals
      INTEGER :: J
      DOUBLE PRECISION :: AP(3),RP(3),E3(3)      
      DO J=1,3 ;AP(J)=P(J)-A(J) ;END DO
      E3(:) = CROSS(E1,AP)
      RP(:) = CROSS(E3,E1)   
      R     = NORM(RP)
      END FUNCTION POINT_TO_LINE_DISTANCE

!---------------------------------------------------------------------
      FUNCTION VECTOR_ANGLE(V1,V2) RESULT(V_ANGLE)
      DOUBLE PRECISION ,INTENT(IN) :: V1(3),V2(3)
      DOUBLE PRECISION :: V_ANGLE
      !Locals
      DOUBLE PRECISION :: H,COS_ANG,SIN_ANG,SMALL ,V3(3)
      SMALL = 1.E-7
      H     = NORM(V1)*NORM(V2)
      IF (H<SMALL) THEN
        V_ANGLE = 0.D0
      ELSE
        COS_ANG = (V1(1)*V2(1)+V1(2)*V2(2)+V1(3)*V2(3))/H
        V3(:)   = CROSS(V1,V2)
        SIN_ANG = NORM(V3)/H
        V_ANGLE = ATAN2(SIN_ANG,COS_ANG) ![RAD]
      END IF      
      END FUNCTION VECTOR_ANGLE

      

!     ******************************************************************
!     *                                                                *
!     *                   SUBROUTINE DGAUSS                            *
!     *                                                                *
!     *         GAUSS ELIMINATION METHOD (SYSTEM SOLUTION)             *
!     *                                                                *
!     ******************************************************************          
   
      SUBROUTINE DGAUSS(ND,N,A,B,IOPT,D,IREG,L)
      
      INTEGER:: IREG ,N ,ND ,NM , IOPT ,K ,KP ,M ,I ,J ,L
      DOUBLE PRECISION:: A,B,D,S,EE
      DIMENSION:: A(ND,ND),B(ND),L(ND) 
       
!     DECOMPOSITION OF THE MATRIX A AND/OR SOLUTION OF THE
!     LINEAR SYSTEM AX=B BY GAUSS' ELIMINATION METHOD

!     INPUT:
!     ND: MAXIMUM COLUMN DIMENSION OF ARRAYS
!     N: DIMENSION OF ARRAYS
!     A: SYSTEM MATRIX (NXN)
!     B: RIGHT-HAND SIDE VECTOR (N)
!     IOPT = 0, DECOMPOSITION OF A AND SOLUTION OF THE SYSTEM
!          = 1, DECOMPOSITION OF A ONLY
!          = 2, SOLUTION OF THE SYSTEM (SKIPS DECOMPOSITION)
        
!     OUTPUT:
!     A: DECOMPOSED MATRIX (LOWER\UPPER PART), IF IOPT=0 OR 1
!     B:=X: SOLUTION VECTOR STORED IN B (N)
!     D: DETERMINANT OF A                                                
!     IREG = 0, SINGULAR MATRIX A
!          = 1, REGULAR MATRIX A
!     L: VECTOR OF PIVOTING INDICES (N)
    
      EE=1.D-12
      IREG=1
      NM=N-1
      IF(IOPT.EQ.2) GO TO 60

!     TRIANGULAR DECOMPOSITION OF A WITH MAXIMUM COLUMN PIVOTING
!     (ELIMINATION WITH MULTIPLIERS STORED IN THE LOWER PART OF A)
	   
      !D=1.D0  !<--- The determinant has great values => Floating Point Overflow error
      DO K=1,NM
         KP=K+1
         M=K
         DO I=KP,N
             IF(DABS(A(I,K)).GT.DABS(A(M,K))) M=I
         End do
         L(K)=M
         !D=D*A(M,K)
         IF(DABS(A(M,K)).LE.EE) GO TO 100   
             IF(M.EQ.K) GO TO 30
                 !D=-D
                 DO J=K,N
                     S=A(M,J)
                     A(M,J)=A(K,J)
                     A(K,J)=S
                 End do
   30        DO I=KP,N
                 S=A(I,K)/A(K,K)
                 A(I,K)=S
                 DO J=KP,N
                     A(I,J)=A(I,J)-S*A(K,J)
                 End do
             End Do
      END Do
      !D=D*A(N,N)
      IF(DABS(A(N,N)).LE.EE) GO TO 100
                                                  
      IF(IOPT.EQ.1) RETURN

!     ELIMINATION OPERATIONS ON B

   60 DO K=1,NM
         KP=K+1    
         M=L(K)
         S=B(M)
         B(M)=B(K)
         B(K)=S
         DO I=KP,N
             B(I)=B(I)-A(I,K)*S
	     END Do
      END DO
      
!     BACK-SUBSTITUTION

      B(N)=B(N)/A(N,N)
      DO I=1,NM
         K=N-I
         KP=K+1
         S=0.D0
         DO J=KP,N
             S=S+A(K,J)*B(J)
         END DO
         B(K)=(B(K)-S)/A(K,K)
      END DO
      RETURN
  100 IREG=0
      RETURN
      
      END SUBROUTINE DGAUSS	 

      END MODULE SINUS_SOLVER 
	  

!     ******************************************************************
!     *                                                                *
!     *                    THE ACTUAL PROGRAM                          *
!     *                                                                *
!     *                                                                *
!     *                                                                *
!     ****************************************************************** 
 
      PROGRAM SOLVER
      USE SINUS_SOLVER
      IMPLICIT NONE
      CALL MAIN
      END PROGRAM SOLVER   


	 


	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  
	 
	 
	 
	 
	 
	 
	 
	 
	 
	 
	 
	 
