<|start_filename|>welib/vortilib/elements/fortran/UIVortexPoint2DPerdiodic.f90<|end_filename|>
!> 
module VortexPoint2DPeriodic
    use PrecisionMod,   only: MK
    implicit none
    real(MK),parameter :: MINNORM2=1e-8
contains
    subroutine fUi_VortexPoint2DPeriodicX_11(DeltaP,Intensity,Lambda,SmoothModel,SmoothParam,bComputeGrad,Uind, Grad)
        use MathConstants , only: twopi
        implicit none 
        ! Input/output arguments 
        real(MK), dimension(3), intent(in) :: DeltaP       !< 3 x 1   Pcp-Pv
        real(MK), intent(in)               :: Intensity    !<  Circulation Gamma
        real(MK), intent(in)               :: Lambda       !<  Period in the X direction
        integer, intent(in) ,optional      :: SmoothModel  !< 
        real(MK), intent(inout),optional      :: SmoothParam  !< 
        logical, intent(in)                :: bComputeGrad !< 
        real(MK), dimension(3),intent(out) :: Uind         !< No Side effects. Note that
        real(MK), dimension(9),intent(out) :: Grad         !< No Side effects
        !         ! Variables declaration 
        real(MK) :: x     !< 
        real(MK) :: y     !< 
        real(MK) :: r2     !< 
        real(MK) :: factor !< 
        real(MK) :: E      !< 
        real(MK) :: rho2   !< 
        real(MK) :: Qm     !< 
        ! 
        ! No need of gradient in 2D
        if(bComputeGrad) then
            Grad(1:9)=0.0_MK !
        else
            Grad(1:9)=0.0_MK
        endif
        factor = 1._MK
        if(SmoothModel>0 .and. SmoothParam>0) then
            rho2   = r2/SmoothParam**2;
            E      = exp(-rho2);
            Qm     = 1._MK;
            factor = (1._MK-Qm*E);
            !
        else
            SmoothParam=0.0_MK
        endif
        Uind=0.0_MK

        r2  = DeltaP(1)**2+ DeltaP(2)**2
        if (r2<MINNORM2) then! to avoid strong singularity
            Uind(3)=0.0_MK
            Grad(1:9)=0.0_MK 
        else
            x=DeltaP(1)*twopi/Lambda
            y=DeltaP(2)*twopi/Lambda
            Uind(1)= -(Intensity*Sinh(y))/(2.*(-Cos(x) + Cosh(y)+SmoothParam**2))
            Uind(2)=  (Intensity* Sin(x))/(2.*(-Cos(x) + Cosh(y)+SmoothParam**2))
            Uind(3)=0.0_MK

        endif

    end subroutine
            
end module VortexPoint2DPeriodic


<|start_filename|>welib/vortilib/elements/fortran/UIPanels.f90<|end_filename|>
!>
module UIPanels
    use PrecisionMod,only:MK
    implicit none
    integer, save:: nWarning=0

contains
    !>
    subroutine fUi_QuadDoubletCst_1n(P1,P2,P3,P4,CPs,Intensity, SmoothModel, SmoothParam,bComputeGrad,bTrailedOnly,UI,Grad,nCPs)
        !use MathConstants,only : pi,fourpi
        use UISegments, only : fUi_SegmentCst_11
        ! Arguments declarations
        integer, intent(in)                        :: nCPs        ! >
        logical, intent(in)                        :: bComputeGrad                     !< 
        logical, intent(in)                        :: bTrailedOnly
        real(MK), dimension(3), intent(in)         :: P1,P2,P3,P4                      !< 3
        real(MK), dimension(3,nCPs), intent(in)    :: CPs                              !< 3 x n
        real(MK), dimension(3,nCPs), intent(inout) :: UI                               !< 3 x n
        real(MK), dimension(9,nCPs), intent(inout) :: Grad                             !< 9 x n
        real(MK), intent(in)                       :: Intensity                        !< 
        integer, intent(in)                        :: SmoothModel        ! >
        real(MK), intent(in)                       :: SmoothParam                  ! >
        ! Variables declarations
        real(MK), dimension(3)                     :: UItmp
        real(MK), dimension(9)                     :: Gradtmp
        real(MK), dimension(3) :: DP1          !< 
        real(MK), dimension(3) :: DP2          !< 
        integer :: i
        ! UItmp, is incremented at each call
        !          write(*,*) 'P1', P1
        !          write(*,*) 'P2', P2
        !          write(*,*) 'P3', P3
        !          write(*,*) 'P4', P4
        do i=1,nCPs
            !              write(*,*) 'CPs', CPs(:,i)
            ! --------------------------------------------------------------------------------
            ! --- "Shed" segments 
            ! --------------------------------------------------------------------------------
            if(.not.bTrailedOnly) then
                ! 1-2 segment
                DP1=CPs(:,i)-P1; DP2=CPs(:,i)-P2;
                call fUi_SegmentCst_11(DP1, DP2,Intensity,SmoothModel,SmoothParam,bComputeGrad,UItmp,Gradtmp)
                UI(:,i)=UI(:,i)+UItmp
                Grad(:,i)=Grad(:,i)+Gradtmp
                !              write(*,*) '1', UItmp, UI(:,i)
                ! 3-4 segment
                DP1=CPs(:,i)-P3; DP2=CPs(:,i)-P4 ;
                call  fUi_SegmentCst_11(DP1,DP2,Intensity,SmoothModel,SmoothParam,bComputeGrad,UItmp,Gradtmp)
                UI(:,i)=UI(:,i)+UItmp
                Grad(:,i)=Grad(:,i)+Gradtmp
                !              write(*,*) '3', UItmp, UI(:,i)
            endif
            ! --------------------------------------------------------------------------------
            ! --- "Trailed" segments 
            ! --------------------------------------------------------------------------------
            ! 2-3 segment
            DP1=CPs(:,i)-P2; DP2=CPs(:,i)-P3;
            call  fUi_SegmentCst_11(DP1, DP2, Intensity, SmoothModel , SmoothParam, bComputeGrad, UItmp,Gradtmp)
            UI(:,i)=UI(:,i)+UItmp
            Grad(:,i)=Grad(:,i)+Gradtmp
            !              write(*,*) '2', UItmp, UI(:,i)
            ! 4-1 segment
            DP1=CPs(:,i)-P4; DP2=CPs(:,i)-P1;
            call  fUi_SegmentCst_11(DP1, DP2, Intensity, SmoothModel , SmoothParam, bComputeGrad, UItmp,Gradtmp)
            UI(:,i)=UI(:,i)+UItmp
            Grad(:,i)=Grad(:,i)+Gradtmp
            !              write(*,*) '4', UItmp, UI(:,i)
            !              UI(:,i)=UItmp
        enddo
    end subroutine  fUi_QuadDoubletCst_1n


    !> 
    subroutine fUi_QuadSourcePlaneCst_11(CP,Sigma,xi,eta,RefPoint,TransfoMat,Viscous,bComputeGrad,UI,Grad)
        use MathConstants,only : pi,fourpi
        use PrecisionMod, only : precision_equal
        real(MK),parameter :: eps_quadsource=1e-6_MK !!!!!!!!!!!!!!!!!! !< Used if z coordinate close to zero
        ! Arguments declarations 
        logical, intent(in)                  :: bComputeGrad        !< 
        real(MK), dimension(3), intent(in)   :: CP                  !< 3 x n
        real(MK), dimension(3), intent(out)  :: UI                  !< 3 x n
        real(MK), dimension(9), intent(out)  :: Grad                !< 9 x n
        real(MK), intent(in)                 :: Sigma              !< 
        real(MK), dimension(3), intent(in)   :: RefPoint            !< 3 x n !coordinate of panel origin in ref coordinates
        real(MK), dimension(4), intent(in)   :: xi!1,xi2,xi3,xi4     !< 4 x n
        real(MK), dimension(4), intent(in)   :: eta!1,eta2,eta3,eta4 !< 4 x n
        real(MK), dimension(3,3), intent(in) :: TransfoMat          !< 3 x 3 x n
        real(MK), intent(in)                 :: Viscous ! TODO
        ! Variable declarations 
        real(MK), dimension(3,3) :: A                          !< 
        real(MK), dimension(3,3) :: tA                         !< 
        real(MK)                 :: d12, d23, d34, d41         !< 
        real(MK)                 :: m12, m23, m34, m41         !< 
        real(MK)                 :: xi1,  xi2,  xi3,  xi4      !< 
        real(MK)                 :: eta1,  eta2,  eta3,  eta4  !< 
        real(MK)                 :: e1,  e2,  e3,  e4          !< 
        real(MK)                 :: h1,  h2,  h3,  h4          !< 
        real(MK)                 :: r1,  r2,  r3,  r4          !< 
        real(MK)                 :: RJ12, RJ23, RJ34, RJ41     !< 
        real(MK)                 :: TAN12, TAN23, TAN34, TAN41 !< 
        real(MK), dimension(3) :: V     !< 
        real(MK), dimension(3) :: P    !< 
        real(MK), dimension(3) :: Pe   !< 
        xi1=xi(1)
        xi2=xi(2)
        xi3=xi(3)
        xi4=xi(4)
        eta1=eta(1)
        eta2=eta(2)
        eta3=eta(3)
        eta4=eta(4)
        ! TODO viscous is not used yet 

        Grad=0.0_MK
        V = 0.0_MK 
        !param that are constant for each panel, distances and slopes - The slopes can be if divided by zero NaN => security required 
        d12=sqrt((xi2-xi1)**2+(eta2-eta1)**2)
        d23=sqrt((xi3-xi2)**2+(eta3-eta2)**2)
        d34=sqrt((xi4-xi3)**2+(eta4-eta3)**2)
        d41=sqrt((xi1-xi4)**2+(eta1-eta4)**2)

        ! transform control points in panel coordinate system using matrix - vectorial way 
        A=TransfoMat(1:3,1:3)
        tA=transpose(A)

        P(1:3)=CP(1:3)-RefPoint(1:3)
        Pe = matmul(A,P) ! transfo in element coordinate system, noted x,y,z, but in fact xi eta zeta
        ! scalars
        r1=sqrt((Pe(1)-xi1)**2 + (Pe(2)-eta1)**2 + Pe(3)**2)
        r2=sqrt((Pe(1)-xi2)**2 + (Pe(2)-eta2)**2 + Pe(3)**2)
        r3=sqrt((Pe(1)-xi3)**2 + (Pe(2)-eta3)**2 + Pe(3)**2)
        r4=sqrt((Pe(1)-xi4)**2 + (Pe(2)-eta4)**2 + Pe(3)**2)

        e1=Pe(3)**2 + (Pe(1)-xi1)**2
        e2=Pe(3)**2 + (Pe(1)-xi2)**2
        e3=Pe(3)**2 + (Pe(1)-xi3)**2
        e4=Pe(3)**2 + (Pe(1)-xi4)**2
        ! 
        h1=(Pe(2)-eta1)*(Pe(1)-xi1)
        h2=(Pe(2)-eta2)*(Pe(1)-xi2)
        h3=(Pe(2)-eta3)*(Pe(1)-xi3)
        h4=(Pe(2)-eta4)*(Pe(1)-xi4)

        ! Velocities in element frame 
        ! --------------------------------------------------------------------------------
        ! --- Log term 
        ! --------------------------------------------------------------------------------
        ! Security - Katz Plotkin page 608 appendix D Code 11 
        if ( r1+r2-d12<=0.0_MK .or. d12 <=0.0_MK ) then
            RJ12=0._MK
        else
            RJ12=1/d12 * log((r1+r2-d12)/(r1+r2+d12))
        endif

        if ( r2+r3-d23<=0.0_MK .or. d23 <=0.0_MK ) then
            RJ23=0._MK
        else
            RJ23=1/d23 * log((r2+r3-d23)/(r2+r3+d23))
        endif

        if ( r3+r4-d34<=0.0_MK .or. d34 <=0.0_MK ) then
            RJ34=0._MK
        else
            RJ34=1/d34 * log((r3+r4-d34)/(r3+r4+d34))
        endif

        if ( r4+r1-d41<=0.0_MK .or. d41 <=0.0_MK ) then
            RJ41=0._MK
        else
            RJ41=1/d41 * log((r4+r1-d41)/(r4+r1+d41))
        endif
        ! --------------------------------------------------------------------------------
        ! --- Tan term 
        ! --------------------------------------------------------------------------------
        ! 12
        if (precision_equal(xi2,xi1)) then ! Security - Hess 1962 - page 47 - bottom
            TAN12=0._MK
        else
            m12=(eta2-eta1)/(xi2-xi1)
            if( abs(Pe(3))<eps_quadsource ) then ! case where z is too small, jumps may occur 
                TAN12=pi*aint((mysign(1.0_MK,(m12*e1-h1)) - mysign(1.0_MK,(m12*e2-h2)))/2) ! Security-Hess1962-page47-top
            else
                TAN12= atan((m12*e1-h1)/(Pe(3)*r1)) - atan((m12*e2-h2)/(Pe(3)*r2))
            endif
        endif
        ! 23
        if (precision_equal(xi3,xi2)) then ! Security - Hess 1962 - page 47 - bottom
            TAN23=0._MK
        else
            m23=(eta3-eta2)/(xi3-xi2)
            if( abs(Pe(3))<eps_quadsource ) then ! case where z is too small, jumps may occur 
                TAN23=pi*aint((mysign(1.0_MK,(m23*e2-h2)) - mysign(1.0_MK,(m23*e3-h3)))/2) ! Security-Hess1962-page47-top
            else
                TAN23= atan((m23*e2-h2)/(Pe(3)*r2)) - atan((m23*e3-h3)/(Pe(3)*r3))
            endif
        endif 
        ! 34
        if (precision_equal(xi4,xi3)) then ! Security - Hess 1962 - page 47 - bottom
            TAN34=0._MK
        else
            m34=(eta4-eta3)/(xi4-xi3)
            if( abs(Pe(3))<eps_quadsource ) then ! case where z is too small, jumps may occur 
                TAN34=pi*aint((mysign(1.0_MK,(m34*e3-h3)) - mysign(1.0_MK,(m34*e4-h4)))/2) ! Security-Hess1962-page47-top
            else
                TAN34= atan((m34*e3-h3)/(Pe(3)*r3)) - atan((m34*e4-h4)/(Pe(3)*r4))
            endif
        endif
        ! 41
        if (precision_equal(xi1,xi4)) then ! Security - Hess 1962 - page 47 - bottom
            TAN41=0._MK
        else
            m41=(eta1-eta4)/(xi1-xi4)
            if( abs(Pe(3))<eps_quadsource ) then ! case where z is too small, jumps may occur 
                TAN41=pi*aint((mysign(1.0_MK,(m41*e4-h4)) - mysign(1.0_MK,(m41*e1-h1)))/2) ! Security-Hess1962-page47-top
            else
                TAN41= atan((m41*e4-h4)/(Pe(3)*r4)) - atan((m41*e1-h1)/(Pe(3)*r1))
            endif
        endif


        ! --------------------------------------------------------------------------------
        ! --- Velocity  in Panel frame
        ! --------------------------------------------------------------------------------
        ! Vx 
        V(1)= Sigma/(fourpi)*( (eta2-eta1)*RJ12 + &
            (eta3-eta2)*RJ23 + &
            (eta4-eta3)*RJ34 + &
            (eta1-eta4)*RJ41 )
        ! vy 

        V(2)= Sigma/(fourpi)*( (xi1-xi2) *RJ12 + &
            (xi2-xi3) *RJ23 + &
            (xi3-xi4) *RJ34 + &
            (xi4-xi1) *RJ41 )

        ! Vz 
        V(3)= Sigma/(fourpi)*( ( TAN12 ) + &
            ( TAN23 ) + &
            ( TAN34 ) + &
            ( TAN41 ) )

        ! --------------------------------------------------------------------------------
        ! --- Velocity in Reference frame 
        ! --------------------------------------------------------------------------------
        UI(1:3)=matmul(tA,V(1:3)) ! summation for all panels
        !              write(*,*) 'A(1,:)',A(1,:)
        !              write(*,*) 'A(2,:)',A(2,:)
        !              write(*,*) 'A(3,:)',A(3,:)
        !              write(*,*) 'TAN12', TAN12
        !              write(*,*) 'TAN23', TAN23
        !              write(*,*) 'TAN34', TAN34
        !              write(*,*) 'TAN41', TAN41
        !              write(*,*) 'RJ', RJ12,RJ23,RJ34,RJ41
        !              write(*,*) 'eta', eta
        !              write(*,*) 'xi', xi
        !              write(*,*) 'V', V
        !              write(*,*) 'UI', UI
        !              if (ip==1) then
        !                  STOP
        !              endif
        !  

        if (bComputeGrad) then
          if (nWarning<2) then
            write(*,*)'Panels: Smooth and Grad TODO'
            nWarning=nWarning+1
          endif
            ! Just a fake test to use Viscous
            if (int(Viscous)==-1000) then
            endif
        endif
    end subroutine  fUi_QuadSourcePlaneCst_11


    !> Same as QuadDoubletCst but uses precomputed length l
    subroutine fUi_QuadDoubletCst_11_b(CP,P1,P2,P3,P4,l2,Intensity,SmoothModel,SmoothParam,bComputeGrad,Uind,Grad)
        use UISegments, only : fUi_SegmentCst_11_b
        ! Arguments declarations 
        real(MK), dimension(3), intent(in)  :: CP           !< 
        real(MK), dimension(3), intent(in)  :: P1           !< 
        real(MK), dimension(3), intent(in)  :: P2           !< 
        real(MK), dimension(3), intent(in)  :: P3           !< 
        real(MK), dimension(3), intent(in)  :: P4           !< 
        real(MK), dimension(4), intent(in)  :: l2           !< Norm2 Length of sides
        real(MK),  intent(in)               :: Intensity    !< 
        integer , intent(in)                :: SmoothModel  !< 
        real(MK) , intent(in)               :: SmoothParam  !< 
        logical, intent(in)                 :: bComputeGrad !< 
        real(MK), dimension(3), intent(out) :: Uind         !< No Side effects!!!
        real(MK), dimension(9), intent(out) :: Grad         !< No Side effects!!!
        ! Variable declarations 
        real(MK), dimension(3) :: Uindtmp            !< No Side effects!!!
        real(MK), dimension(9) :: Gradtmp            !< No Side effects!!!
        real(MK), dimension(3) :: DP1, DP2
        ! Side1 - 1-2 segment
        DP1=CP-P1; DP2=CP-P2
        call fUi_SegmentCst_11_b( DP1, DP2, l2(1), Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uindtmp(1:3) ! Important first pass
        Grad(1:9)=Gradtmp(1:9)
        !print "(A,3F21.16)","ui :",Uindtmp
        ! Side2 - 2-3 segment
        DP1=CP-P2; DP2=CP-P3
        call fUi_SegmentCst_11_b( DP1 , DP2, l2(2), Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uind(1:3)+Uindtmp(1:3)
        Grad(1:9)=Grad(1:9)+Gradtmp(1:9)
        !print "(A,3F21.16)","ui :",Uindtmp
        ! Side3 -3-4 segment
        DP1=CP-P3; DP2=CP-P4
        call fUi_SegmentCst_11_b( DP1 , DP2, l2(3), Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uind(1:3)+Uindtmp(1:3)
        Grad(1:9)=Grad(1:9)+Gradtmp(1:9)
        !print "(A,3F21.16)","ui :",Uindtmp
        ! Side4 -4-1 segment
        DP1=CP-P4; DP2=CP-P1
        call fUi_SegmentCst_11_b( DP1, DP2, l2(4), Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uind(1:3)+Uindtmp(1:3)
        Grad(1:9)=Grad(1:9)+Gradtmp(1:9)
        !print "(A,3F21.16)","ui :",Uindtmp
    end subroutine  fUi_QuadDoubletCst_11_b

    !> Same as HorseShoe_11 but uses precomputed length l
    subroutine fUi_HorseShoe_11_b(CP,P1,P2,P3,P4,l2,Intensity,SmoothModel,SmoothParam,bComputeGrad,Uind,Grad)
        use UISegments, only : fUi_SegmentCst_11_b
        ! Arguments declarations 
        real(MK), dimension(3), intent(in)  :: CP           !< 
        real(MK), dimension(3), intent(in)  :: P1           !< 
        real(MK), dimension(3), intent(in)  :: P2           !< 
        real(MK), dimension(3), intent(in)  :: P3           !< 
        real(MK), dimension(3), intent(in)  :: P4           !< 
        real(MK), dimension(4), intent(in)  :: l2           !< Norm2 Length of sides
        real(MK),  intent(in)               :: Intensity    !< 
        integer , intent(in)                :: SmoothModel  !< 
        real(MK) , intent(in)               :: SmoothParam  !< 
        logical, intent(in)                 :: bComputeGrad !< 
        real(MK), dimension(3), intent(out) :: Uind         !< No Side effects!!!
        real(MK), dimension(9), intent(out) :: Grad         !< No Side effects!!!
        ! Variable declarations 
        real(MK), dimension(3) :: Uindtmp            !< No Side effects!!!
        real(MK), dimension(9) :: Gradtmp            !< No Side effects!!!
        real(MK), dimension(3) :: DP1          !< 
        real(MK), dimension(3) :: DP2          !< 
        ! Side1 - 1-2 segment
        DP1=CP-P1 ; DP2= CP-P2;
        call fUi_SegmentCst_11_b( DP1, DP2, l2(1), Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uindtmp(1:3) ! Important first pass
        Grad(1:9)=Gradtmp(1:9)
        !print "(A,3F21.16)","ui :",Uindtmp
        ! Side2 - 2-3 segment
        DP1=CP-P2 ; DP2= CP-P3;
        call fUi_SegmentCst_11_b( DP1 , DP2, l2(2), Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uind(1:3)+Uindtmp(1:3)
        Grad(1:9)=Grad(1:9)+Gradtmp(1:9)
        !print "(A,3F21.16)","ui :",Uindtmp
        ! Side4 -4-1 segment
        DP1=CP-P4 ; DP2= CP-P1;
        call fUi_SegmentCst_11_b( DP1, DP2, l2(4), Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uind(1:3)+Uindtmp(1:3)
        Grad(1:9)=Grad(1:9)+Gradtmp(1:9)
        !print "(A,3F21.16)","ui :",Uindtmp
    end subroutine  fUi_HorseShoe_11_b
    !>  
    subroutine fUi_QuadDoubletCst_11(CP,P1,P2,P3,P4,Intensity,SmoothModel,SmoothParam,bComputeGrad,Uind,Grad)
        use UISegments, only : fUi_SegmentCst_11
        ! Arguments declarations 
        real(MK), dimension(3), intent(in)  :: CP           !< 
        real(MK), dimension(3), intent(in)  :: P1           !< 
        real(MK), dimension(3), intent(in)  :: P2           !< 
        real(MK), dimension(3), intent(in)  :: P3           !< 
        real(MK), dimension(3), intent(in)  :: P4           !< 
        real(MK),  intent(in)               :: Intensity    !< 
        integer , intent(in)                :: SmoothModel  !< 
        real(MK) , intent(in)               :: SmoothParam  !< 
        logical, intent(in)                 :: bComputeGrad !< 
        real(MK), dimension(3), intent(out) :: Uind         !< No Side effects!!!
        real(MK), dimension(9), intent(out) :: Grad         !< No Side effects!!!
        ! Variable declarations 
        real(MK), dimension(3) :: Uindtmp    !< No Side effects!!!
        real(MK), dimension(9) :: Gradtmp    !< No Side effects!!!
        real(MK), dimension(3) :: DP1, DP2
        ! Side1 - 1-2 segment
        DP1=CP-P1; DP2=CP-P2
        call fUi_SegmentCst_11 ( DP1, DP2, Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uindtmp(1:3) ! Important first pass
        Grad(1:9)=Gradtmp(1:9)
!         print "(A,3F35.30)","ui :",Uindtmp
        ! Side2 - 2-3 segment
        DP1=CP-P2; DP2=CP-P3
        call fUi_SegmentCst_11 ( DP1, DP2, Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uind(1:3)+Uindtmp(1:3)
        Grad(1:9)=Grad(1:9)+Gradtmp(1:9)
!         print "(A,3F35.30)","ui :",Uindtmp
        ! Side3 - 3-4 segment
        DP1=CP-P3; DP2=CP-P4
        call fUi_SegmentCst_11 ( DP1, DP2, Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uind(1:3)+Uindtmp(1:3)
        Grad(1:9)=Grad(1:9)+Gradtmp(1:9)
!         print "(A,3F35.30)","ui :",Uindtmp
        ! Side4 - 4-1 segment
        DP1=CP-P4; DP2=CP-P1
        call fUi_SegmentCst_11 ( DP1, DP2, Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uind(1:3)+Uindtmp(1:3)
        Grad(1:9)=Grad(1:9)+Gradtmp(1:9)
!         print "(A,3F35.30)","ui :",Uindtmp
    end subroutine  fUi_QuadDoubletCst_11


    !> Like constantdoublet /vortex ring but without the 34 segment
    subroutine fUi_HorseShoe_11(CP,P1,P2,P3,P4,Intensity,SmoothModel,SmoothParam,bComputeGrad,Uind,Grad)
        use UISegments, only : fUi_SegmentCst_11
        ! Arguments declarations 
        real(MK), dimension(3), intent(in)  :: CP           !< 
        real(MK), dimension(3), intent(in)  :: P1           !< 
        real(MK), dimension(3), intent(in)  :: P2           !< 
        real(MK), dimension(3), intent(in)  :: P3           !< 
        real(MK), dimension(3), intent(in)  :: P4           !< 
        real(MK),  intent(in)               :: Intensity    !< 
        integer , intent(in)                :: SmoothModel  !< 
        real(MK) , intent(in)               :: SmoothParam  !< 
        logical, intent(in)                 :: bComputeGrad !< 
        real(MK), dimension(3), intent(out) :: Uind         !< No Side effects!!!
        real(MK), dimension(9), intent(out) :: Grad         !< No Side effects!!!
        ! Variable declarations 
        real(MK), dimension(3) :: Uindtmp    !< No Side effects!!!
        real(MK), dimension(9) :: Gradtmp    !< No Side effects!!!
        real(MK), dimension(3) :: DP1, DP2
        ! Side1 - 1-2 segment
        DP1=CP-P1; DP2=CP-P2
        call fUi_SegmentCst_11 ( DP1, DP2, Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uindtmp(1:3) ! Important first pass
        Grad(1:9)=Gradtmp(1:9)
!         print "(A,3F35.30)","ui :",Uindtmp
        ! Side2 - 2-3 segment
        DP1=CP-P2; DP2=CP-P3
        call fUi_SegmentCst_11 ( DP1, DP2, Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uind(1:3)+Uindtmp(1:3)
        Grad(1:9)=Grad(1:9)+Gradtmp(1:9)
!         print "(A,3F35.30)","ui :",Uindtmp
        ! Side4 - 4-1 segment
        DP1=CP-P4; DP2=CP-P1
        call fUi_SegmentCst_11 ( DP1, DP2, Intensity , SmoothModel , SmoothParam , bComputeGrad , Uindtmp , Gradtmp ) 
        Uind(1:3)=Uind(1:3)+Uindtmp(1:3)
        Grad(1:9)=Grad(1:9)+Gradtmp(1:9)
!         print "(A,3F35.30)","ui :",Uindtmp
    end subroutine  fUi_HorseShoe_11



    !> 
    subroutine fUi_QuadSourcePlaneCst(CPs,Sigmas,xi,eta,RefPoint,TransfoMat,Viscous,bComputeGrad,UI,Grad,nCPs,nPanels)
        use MathConstants,only : pi,fourpi
        use PrecisionMod, only : precision_equal
        real(MK),parameter :: eps_quadsource=1e-6_MK !!!!!!!!!!!!!!!!!! !< Used if z coordinate close to zero
        ! Arguments declarations 
        logical, intent(in)                          :: bComputeGrad !< 
        integer, intent(in)                          :: nCPs
        integer, intent(in)                          :: nPanels
        real(MK), dimension(3,nCPs), intent(in)      :: CPs            !< 3 x n
        real(MK), dimension(3,nCPs), intent(inout)   :: UI             !< 3 x n
        real(MK), dimension(9,nCPs), intent(inout)   :: Grad           !< 9 x n
        real(MK), dimension(nPanels), intent(in)     :: Sigmas         !< 
        real(MK), dimension(3,nPanels), intent(in)   :: RefPoint       !< 3 x n !coordinate of panel origin in ref coordinates
        real(MK), dimension(4,nPanels), intent(in)   :: xi             !< 4 x n
        real(MK), dimension(4,nPanels), intent(in)   :: eta            !< 4 x n
        real(MK), dimension(3,3,nPanels), intent(in) :: TransfoMat     !< 3 x 3 x n
        real(MK), dimension(:), intent(in)           :: Viscous ! TODO
        ! Variable declarations 
        real(MK), dimension(3,3) :: A  !<  
        real(MK), dimension(3,3) :: tA  !<  
        real(MK) :: d12  !<  
        real(MK) :: d23  !<  
        real(MK) :: d34  !<  
        real(MK) :: d41  !<  
        real(MK) :: m12  !<  
        real(MK) :: m23  !<  
        real(MK) :: m34  !<  
        real(MK) :: m41  !<  
        real(MK) :: e1    !< 
        real(MK) :: e2    !< 
        real(MK) :: e3    !< 
        real(MK) :: e4    !< 
        real(MK) :: h1    !< 
        real(MK) :: h2    !< 
        real(MK) :: h3    !< 
        real(MK) :: h4    !< 
        real(MK) :: r1    !< 
        real(MK) :: r2    !< 
        real(MK) :: r3    !< 
        real(MK) :: r4    !< 
        real(MK) :: RJ12  !< 
        real(MK) :: RJ23  !< 
        real(MK) :: RJ34  !< 
        real(MK) :: RJ41  !< 
        real(MK) :: TAN12 !< 
        real(MK) :: TAN23 !< 
        real(MK) :: TAN34 !< 
        real(MK) :: TAN41 !< 
        real(MK), dimension(3) :: V     !< 
        real(MK), dimension(3) :: P    !< 
        real(MK), dimension(3) :: Pe   !< 
        ! 
        integer :: ip,icp ! loop index

        ! [ ui vi wi] = fUi_QuadSourcePlaneConstantN(1:4,4:7,7:10,0:5,0:5,0:5,0:5,0:5,0:5,0) 


        ! TODO viscous is not used yet 

        V = 0.0_MK 
        do ip=1,nPanels !loop on panels 
            !param that are constant for each panel, distances and slopes - The slopes can be if divided by zero NaN => security required 
            d12=sqrt((xi(2,ip)-xi(1,ip))**2+(eta(2,ip)-eta(1,ip))**2)
            d23=sqrt((xi(3,ip)-xi(2,ip))**2+(eta(3,ip)-eta(2,ip))**2)
            d34=sqrt((xi(4,ip)-xi(3,ip))**2+(eta(4,ip)-eta(3,ip))**2)
            d41=sqrt((xi(1,ip)-xi(4,ip))**2+(eta(1,ip)-eta(4,ip))**2)

            do icp=1,nCPs ! loop on Control Points
                ! transform control points in panel coordinate system using matrix - vectorial way 
                A=TransfoMat(1:3,1:3,ip)
                tA=transpose(A)

                P(1:3)=CPs(1:3,icp)-RefPoint(1:3,ip)
                !             write(*,*) 'Px', P(1,:)
                !             write(*,*) 'Py', P(2,:)
                !             write(*,*) 'Pz', P(3,:)
                Pe = matmul(A,P) ! transfo in element coordinate system, noted x,y,z, but in fact xi eta zeta
                !             write(*,*) 'Pex', Pe(1,:)
                !             write(*,*) 'Pey', Pe(2,:)
                !             write(*,*) 'Pez', Pe(3,:)
                !             write(*,*) 'P', shape(MP)
                !             write(*,*) 'Pe', shape(MPe)
                !             write(*,*) 'r1', shape(r1)
                !             write(*,*) 'expr', shape(sqrt((MPe(1,:)-xi(1,ip))**2))

                ! scalars
                r1=sqrt((Pe(1)-xi(1,ip))**2 + (Pe(2)-eta(1,ip))**2 + Pe(3)**2)
                r2=sqrt((Pe(1)-xi(2,ip))**2 + (Pe(2)-eta(2,ip))**2 + Pe(3)**2)
                r3=sqrt((Pe(1)-xi(3,ip))**2 + (Pe(2)-eta(3,ip))**2 + Pe(3)**2)
                r4=sqrt((Pe(1)-xi(4,ip))**2 + (Pe(2)-eta(4,ip))**2 + Pe(3)**2)

                e1=Pe(3)**2 + (Pe(1)-xi(1,ip))**2
                e2=Pe(3)**2 + (Pe(1)-xi(2,ip))**2
                e3=Pe(3)**2 + (Pe(1)-xi(3,ip))**2
                e4=Pe(3)**2 + (Pe(1)-xi(4,ip))**2
                ! 
                h1=(Pe(2)-eta(1,ip))*(Pe(1)-xi(1,ip))
                h2=(Pe(2)-eta(2,ip))*(Pe(1)-xi(2,ip))
                h3=(Pe(2)-eta(3,ip))*(Pe(1)-xi(3,ip))
                h4=(Pe(2)-eta(4,ip))*(Pe(1)-xi(4,ip))

                ! Velocities in element frame 
                !             write(*,*) 'RJbefore', RJ12,RJ23,RJ34,RJ41
                ! --------------------------------------------------------------------------------
                ! --- Log term 
                ! --------------------------------------------------------------------------------
                ! Security - Katz Plotkin page 608 appendix D Code 11 
                if ( r1+r2-d12<=0.0_MK .or. d12 <=0.0_MK ) then
                    RJ12=0._MK
                else
                    RJ12=1/d12 * log((r1+r2-d12)/(r1+r2+d12))
                endif

                if ( r2+r3-d23<=0.0_MK .or. d23 <=0.0_MK ) then
                    RJ23=0._MK
                else
                    RJ23=1/d23 * log((r2+r3-d23)/(r2+r3+d23))
                endif

                if ( r3+r4-d34<=0.0_MK .or. d34 <=0.0_MK ) then
                    RJ34=0._MK
                else
                    RJ34=1/d34 * log((r3+r4-d34)/(r3+r4+d34))
                endif

                if ( r4+r1-d41<=0.0_MK .or. d41 <=0.0_MK ) then
                    RJ41=0._MK
                else
                    RJ41=1/d41 * log((r4+r1-d41)/(r4+r1+d41))
                endif
                !             write(*,*) 'RJafter', RJ12,RJ23,RJ34,RJ41
                ! 
                ! --------------------------------------------------------------------------------
                ! --- Tan term 
                ! --------------------------------------------------------------------------------
                ! 12
                if (precision_equal(xi(2,ip),xi(1,ip))) then ! Security - Hess 1962 - page 47 - bottom
                    TAN12=0._MK
                else
                    m12=(eta(2,ip)-eta(1,ip))/(xi(2,ip)-xi(1,ip))
                    if( abs(Pe(3))<eps_quadsource ) then ! case where z is too small, jumps may occur 
                        TAN12=pi*aint((mysign(1.0_MK,(m12*e1-h1)) - mysign(1.0_MK,(m12*e2-h2)))/2) ! Security-Hess1962-page47-top
                    else
                        TAN12= atan((m12*e1-h1)/(Pe(3)*r1)) - atan((m12*e2-h2)/(Pe(3)*r2))
                    endif
                endif
                ! 23
                if (precision_equal(xi(3,ip),xi(2,ip))) then ! Security - Hess 1962 - page 47 - bottom
                    TAN23=0._MK
                else
                    m23=(eta(3,ip)-eta(2,ip))/(xi(3,ip)-xi(2,ip))
                    if( abs(Pe(3))<eps_quadsource ) then ! case where z is too small, jumps may occur 
                        TAN23=pi*aint((mysign(1.0_MK,(m23*e2-h2)) - mysign(1.0_MK,(m23*e3-h3)))/2) ! Security-Hess1962-page47-top
                    else
                        TAN23= atan((m23*e2-h2)/(Pe(3)*r2)) - atan((m23*e3-h3)/(Pe(3)*r3))
                    endif
                endif 
                ! 34
                if (precision_equal(xi(4,ip),xi(3,ip))) then ! Security - Hess 1962 - page 47 - bottom
                    TAN34=0._MK
                else
                    m34=(eta(4,ip)-eta(3,ip))/(xi(4,ip)-xi(3,ip))
                    if( abs(Pe(3))<eps_quadsource ) then ! case where z is too small, jumps may occur 
                        TAN34=pi*aint((mysign(1.0_MK,(m34*e3-h3)) - mysign(1.0_MK,(m34*e4-h4)))/2) ! Security-Hess1962-page47-top
                    else
                        TAN34= atan((m34*e3-h3)/(Pe(3)*r3)) - atan((m34*e4-h4)/(Pe(3)*r4))
                    endif
                endif
                ! 41
                if (precision_equal(xi(1,ip),xi(4,ip))) then ! Security - Hess 1962 - page 47 - bottom
                    TAN41=0._MK
                else
                    m41=(eta(1,ip)-eta(4,ip))/(xi(1,ip)-xi(4,ip))
                    if( abs(Pe(3))<eps_quadsource ) then ! case where z is too small, jumps may occur 
                        TAN41=pi*aint((mysign(1.0_MK,(m41*e4-h4)) - mysign(1.0_MK,(m41*e1-h1)))/2) ! Security-Hess1962-page47-top
                    else
                        TAN41= atan((m41*e4-h4)/(Pe(3)*r4)) - atan((m41*e1-h1)/(Pe(3)*r1))
                    endif
                endif


                ! --------------------------------------------------------------------------------
                ! --- Velocity  in Panel frame
                ! --------------------------------------------------------------------------------
                ! Vx 
                V(1)= Sigmas(ip)/(fourpi)*( (eta(2,ip)-eta(1,ip))*RJ12 + &
                    (eta(3,ip)-eta(2,ip))*RJ23 + &
                    (eta(4,ip)-eta(3,ip))*RJ34 + &
                    (eta(1,ip)-eta(4,ip))*RJ41 )

                !          write(*,*)                     (eta(2,ip)-eta(1,ip))*RJ12 
                !          write(*,*)                     (eta(3,ip)-eta(2,ip))*RJ23 
                !          write(*,*)                     (eta(4,ip)-eta(3,ip))*RJ34 
                !          write(*,*)                     (eta(1,ip)-eta(4,ip))*RJ41 
                !                  write(*,*) Sigmas(ip)/(fourpi) 
                !                  write(*,*) Sigmas(ip)
                !                  write(*,*) fourpi

                ! vy 

                V(2)= Sigmas(ip)/(fourpi)*( (xi(1,ip)-xi(2,ip)) *RJ12 + &
                    (xi(2,ip)-xi(3,ip)) *RJ23 + &
                    (xi(3,ip)-xi(4,ip)) *RJ34 + &
                    (xi(4,ip)-xi(1,ip)) *RJ41 )

                ! Vz 
                V(3)= Sigmas(ip)/(fourpi)*( ( TAN12 ) + &
                    ( TAN23 ) + &
                    ( TAN34 ) + &
                    ( TAN41 ) )

                ! --------------------------------------------------------------------------------
                ! --- Velocity in Reference frame 
                ! --------------------------------------------------------------------------------
                UI(1:3,icp)=UI(1:3,icp)+matmul(tA,V(1:3)) ! summation for all panels
                !              write(*,*) 'A(1,:)',A(1,:)
                !              write(*,*) 'A(2,:)',A(2,:)
                !              write(*,*) 'A(3,:)',A(3,:)
                !              write(*,*) 'TAN12', TAN12
                !              write(*,*) 'TAN23', TAN23
                !              write(*,*) 'TAN34', TAN34
                !              write(*,*) 'TAN41', TAN41
                !              write(*,*) 'RJ', RJ12,RJ23,RJ34,RJ41
                !              write(*,*) 'eta', eta
                !              write(*,*) 'xi', xi
                !              write(*,*) 'V', V
                !              write(*,*) 'UI', UI
                !              if (ip==1) then
                !                  STOP
                !              endif
                !  

            end do ! control points
        end do !end panel loop
        if (bComputeGrad) then
            Grad=Grad+0.0_MK
            write(*,*)'Panels: Smooth and Grad TODO'
            ! Just a fake test to use Viscous
            if (int(Viscous(1))==-1000) then
            endif
        endif
        !         write(*,*) UI
        !         ui=Ui(:,transpose(1))
        !         vi=Ui(:,transpose(2))
        !         wi=Ui(:,transpose(3))
    end subroutine  fUi_QuadSourcePlaneCst 






    elemental real(MK) function mysign(ref,val)
        use PrecisionMod, only: PRECISION_EPS
        real(MK),intent(in) ::ref
        real(MK),intent(in) ::val
        if ( abs(val)>PRECISION_EPS ) then
            mysign=sign(ref,val)
        else
            mysign=1.0_MK
        endif
    endfunction

end module UIPanels




<|start_filename|>welib/airfoils/matlab/KarmanTrefftzGeometry.f90<|end_filename|>
program mapsl
    !  Based on: <NAME> December, 2000
    !------------------------------------------------------------------------------­

    !	Declare the variables----------------------------------------------------­
    implicit none
    character(len=1) yesno
    character(len=20) :: buffer
    character(len=20) :: params
    real, parameter :: pi=3.1415927e00, rad=1.7453293e-02
    real :: x,y,theta,rsq,u,v,r,xc,yc,alpha_deg,alpha,rc,beta,rcsq,g,gamma, &
        lambda,tau_deg,xi,eta,x_circle,y_circle,fl2,u_foil,v_foil, &
        theta1,theta2,dtheta,rmax,dr,cp,xle,xte,chord,cl
    real, allocatable, dimension(:) :: xi_foil,cp_foil,eta_foil
    complex :: z,zp1,zm1,zeta,w_circle,dzeta_dz,w_foil
    integer :: n,m,ntheta,nradial

    !input the variables defining foil geometry and angle of attack ----------­
    write(*,'(a)') ' enter xc,yc,trailing edge angle ' ! position of center of circle--------­
    read(*,*) xc,yc,tau_deg ! note: xc must be <= 0 ------------­
    lambda=2.0-tau_deg/180.0 ! the exponent in the k-t mapping function -­
    fl2=4.0*lambda**2 ! pre-comute a mapping function constant----­
    rcsq=(1.0-xc)**2+yc**2 ! radius of circle passing through (1,0) ---­

    rc=sqrt(rcsq) 
    beta=atan(yc/(1.0-xc)) 

    alpha_deg=0
    theta1=0
    theta2=360
    dtheta=1
    ntheta=nint((theta2+theta1)/dtheta)+1 ! number of radial grid lines 
    !	alocate the arrays and open output files --------------------------------­
    allocate(xi_foil(ntheta),eta_foil(ntheta),cp_foil(ntheta)) ! arrays for foil cp -----­

    buffer=leading_zero(xc)
    params='_'//trim(buffer)
    buffer=leading_zero(yc)
    params=trim(params)//'_'//trim(buffer)
    write(buffer,'(I0)')int(tau_deg)
    params=trim(params)//'_'//trim(buffer)
!     write(buffer,'(I0)')int(alpha_deg)
!     params=trim(params)//'_'//trim(buffer)

!     open(2,file='cp-karman-trefftz'//trim(params)//'.dat',status='unknown',form='formatted')
    open(3,file='geom-karman-trefftz'//trim(params)//'.dat',status='unknown',form='formatted')

    !-----generate velocity and presure field -------------------------------------­
    m=1
    r=rc+real(m-1)*dr
    rsq=r**2
    ! grid radial coordinate in the z (circle) plane--­
    do n=1,ntheta
        theta=rad*(theta1-real(n-1)*dtheta) ! grid angular coordinate--­
        x=xc+r*cos(theta) ! convert to cartesian-----­
        y=yc+r*sin(theta)
        !-----------compute the velocity field around the circle-----------------------­
        u=cos(alpha)-rcsq*cos(2.0*theta-alpha)/rsq-g*sin(theta)/r
        v=sin(alpha)-rcsq*sin(2.0*theta-alpha)/rsq+g*cos(theta)/r

        !-----------express the field point position and velocity in complex form------­
        z=cmplx(x,y) ! form complex number (x+iy)-------­
        w_circle=cmplx(u,-v)
        ! complex velocity in z plane------­
        !-----------use the karman-trefftz transformation to map points to zeta plane--­
        zp1=(z+1.0)**lambda
        zm1=(z-1.0)**lambda
        zeta=lambda*(zp1+zm1)/(zp1-zm1)
        xi=real(zeta) ! transformed x coordinate---------­
        eta=aimag(zeta) ! transformed y coordinate --------­
        !-----------compute the derivative of the mapping function; transform velocities
        dzeta_dz=fl2*((z-1.0)**(lambda-1.0)*(z+1.0)**(lambda-1.0))/ (zp1-zm1)**2

        if(abs(dzeta_dz)>0.0001) then
            w_foil=w_circle/dzeta_dz ! complex velocity in foil plane------­

        else
            w_foil=0.0
        end if
        u_foil=real(w_foil)
        v_foil=-aimag(w_foil)

        ! avoids divide by zero at trailing edge

        !-----------compute the presure coefficient, cp, and output to plotting file---­

        cp=abs(w_foil)**2-1.0

        if(m==1) then
            ! save presure on foil surface for later use--­
            xi_foil(n)=xi
            eta_foil(n)=eta
            cp_foil(n)=cp
        end if
    end do
    !-----scale the chordwise coordinate to (0,1) and output foil pressure dist.---­
    xle=10.0
    xte=-10.0
    do n=1,ntheta
        xle=min(xle,xi_foil(n))
        xte=max(xte,xi_foil(n))
    end do
    ! Coordinate scalling
    chord=xte-xle
    do n=1,ntheta
        xi_foil(n)=(xi_foil(n)-xle)/chord
        eta_foil(n)=eta_foil(n)/chord
    end do
!     write(2,'(2f10.5)') (xi_foil(n),cp_foil(n),n=1,ntheta) 
!     close(2) 

    write(3,'(2f10.5)') (xi_foil(n),eta_foil(n),n=1,ntheta) 
    close(3) 

contains

    function leading_zero(x) result(buffer)
        real :: x
        character(len=20) :: buffer
        integer :: ix
        if (x.lt.0) then
            ix = ceiling(x)
        else
            ix = floor(x)
        endif
        if (x.lt.0) then
            write(buffer,1) abs(ix), abs(x)-floor(abs(x))
            1   format('-',i0,f0.1)    
        else
            write(buffer,2) floor(x), abs(x)-floor(abs(x))
            2   format(i0,f0.1)
        endif

    end function




end program mapsl


<|start_filename|>welib/vortilib/elements/fortran/PrecisionMod.f90<|end_filename|>
module PrecisionMod
    use iso_c_binding!OTHER-COMPILER
    implicit none
    integer, parameter :: MK = C_DOUBLE  !< MK stand for MyKind : The kind used throughout the program.
    ! Iso c binding pseudo replacement... compiler and platform dependent!
    ! Gfortran and intel have those values, but better to comment and use iso_c_binding above)
    !integer, parameter :: C_INT    = 4 !COMPAQ-COMPILER
    !integer, parameter :: C_FLOAT  = 4 !COMPAQ-COMPILER
    !integer, parameter :: C_DOUBLE = 8 !COMPAQ-COMPILER
    !integer, parameter :: C_CHAR   = 1 !COMPAQ-COMPILER

    integer, parameter :: SP = kind( 1.0e0 ) !< Single Precision
    integer, parameter :: DP = kind( 1.0d0 ) !< Double Precision
    !integer, parameter :: DP = selected_real_kind( 16 ) !< Double Precision
!     integer, parameter :: QP = selected_real_kind (32) !< Quad Precision
    integer, parameter :: QP = selected_real_kind(33, 4931)
    !integer, parameter :: MK = C_DOUBLE       !< MK stand for MyKind : The kind used throughout the program. MODIFIED BY MAKEFILE
    integer, parameter :: CDK = C_DOUBLE       !< C Double Kind
    integer, parameter :: CIK = C_INT          !< C Integer Kind
    integer, parameter :: CIC = C_CHAR         !< C Char Kind

    real(MK),parameter :: PRECISION_EPS =  epsilon(1.0_MK) !< Machine Precision For the given MyKind for problems of scale 1!
!     real(MK),parameter :: PRECISION_UI  = 2.2204460e-10_MK  !< 
    real(MK),parameter :: PRECISION_UI  = PRECISION_EPS/100  !< 
    ! #define FLT_EPSILON 1.19209290e-07F
    ! #define DBL_EPSILON 2.2204460492503131e-16
    ! #define LDBL_EPSILON 1.0842021724855044340075E-19L
    ! #define PRECISION_EPS DBL_EPSILON //< Machine Precision For the given MyKind for problems of scale 1!
   
    ! Pointers
    !integer, parameter :: IPTR_KIND = int_ptr_kind() !INTEL-COMPILER !For compatibility 32/64bits 
    integer, parameter :: IPTR_KIND = 4               !OTHER-COMPILER !For compatibility 32/64bits 
    !integer, parameter :: IPTR_KIND = 4              !COMPAQ-COMPILER!For compatibility 32/64bits


    integer, parameter :: I4 = 4
    
contains

    !> Checks if two reals are equal given machine precision
    logical function precision_equal(x,y)
        real(MK), intent(in) :: x
        real(MK), intent(in) :: y
        precision_equal=abs(x-y)<PRECISION_EPS
    end function

    !> Checks if two reals are different
    logical function precision_different(x,y)
        real(MK), intent(in) :: x
        real(MK), intent(in) :: y
        precision_different=.not.precision_equal(x,y)
    end function
    
    !> Checks if two reals are equal given machine precision
    logical function precision_equal_soft(x,y)
        real(MK), intent(in) :: x
        real(MK), intent(in) :: y
        precision_equal_soft=abs(x-y)<1000*PRECISION_EPS
    end function

    !> Checks if two reals are equal given machine precision
    logical function precision_equal_r4(x,y)
        real(4), intent(in) :: x
        real(4), intent(in) :: y
        precision_equal_r4=abs(x-y)<epsilon(real(1.0,4))
    end function

    !> Checks if two reals are different
    logical function precision_different_r4(x,y)
        real(4), intent(in) :: x
        real(4), intent(in) :: y
        precision_different_r4=.not.precision_equal_r4(x,y)
    end function

    !> Checks if two reals are equal given machine precision
    logical function precision_equal_r8(x,y)
        real(8), intent(in) :: x
        real(8), intent(in) :: y
        precision_equal_r8=abs(x-y)<epsilon(real(1.0,8))
    end function

    !> Checks if two reals are different
    logical function precision_different_r8(x,y)
        real(8), intent(in) :: x
        real(8), intent(in) :: y
        precision_different_r8=.not.precision_equal_r8(x,y)
    end function

    !> Checks if two reals are equal given machine precision
    logical function precision_equal_rdp(x,y)
        double precision, intent(in) :: x
        double precision, intent(in) :: y
        precision_equal_rdp=abs(x-y)<epsilon(1.0d0)
    end function

    !> Checks if two reals are different
    logical function precision_different_rdp(x,y)
        double precision, intent(in) :: x
        double precision, intent(in) :: y
        precision_different_rdp=.not. precision_equal_rdp(x,y)
    end function

    !> Print precision 
    subroutine precision_print()
        print'(A,A)'     ,'Precision : ',precision_string()
        print'(A,5I3)'   ,'Precision types: ',MK,SP,DP,C_FLOAT,C_DOUBLE
        print'(a,E80.64)','Precision eps: ',PRECISION_EPS
        print'(a,E80.64)','Precision ui : ',PRECISION_UI
        print'(A,I3)','Bytes MK: ',sizeof(1.0_MK)!OTHER-COMPILER
        print'(A,I3)','Bytes SP: ',sizeof(1.0_SP)!OTHER-COMPILER
        print'(A,I3)','Bytes DP: ',sizeof(1.0_DP)!OTHER-COMPILER
    end subroutine
    
    integer function sizeOfReal() result(num_bytes)
        !num_bytes=0             !COMPAQ-COMPILER
        num_bytes=int(sizeof(1.0_MK))!OTHER-COMPILER
!         implicit none
!         character, dimension(:), allocatable :: c
!         real(MK) :: r
!         r=0.0_MK
! !         intrinsic size
! !         intrinsic transfer
! !         sizeOfReal = size(transfer(r,c))
!         inquire(iolength=sizeOfReal)r
!          use iso_fortran_env
!         integer :: num_file_storage_units
!         inquire(iolength=num_file_storage_units) 1.0D0
!         num_bytes = num_file_storage_units*FILE_STORAGE_SIZE/8
!         write(*,*) "double has size: ", num_bytes
!         return
    end function sizeOfReal 


    !> Print precision 
    function precision_string() result(s)
        character(len=40) :: s
        if(MK==SP) then
            if(SP==C_FLOAT) then
                s="single precision - C float"
            else
                s="single precision"
            endif
        elseif(MK==DP) then
            if(SP==C_DOUBLE) then
                s="double precision - C double"
            else
                s="double precision"
            endif
        else
            if(MK==C_FLOAT) then
            s="C float"
            elseif(MK==C_DOUBLE) then
                s="C double"
            else
                s="unknown"
            endif
        endif
    end function

end module PrecisionMod


<|start_filename|>welib/vortilib/elements/fortran/UISegments.f90<|end_filename|>
!> Segments Induced velocity functions.
module UISegments
    use PrecisionMod, only: MK,PRECISION_UI
    implicit none

    ! Norm for the simple test exactly on singulartiy
    real(MK),parameter :: MINNORMSIMP=1e-6_MK
    ! Min norm for escaping if denominator or cross prod too small ! Has quite some influence!!!
    ! real(MK),parameter :: MINNORM=1e-08_MK
    ! real(MK),parameter :: MINNORM2=1e-16_MK
    ! Added to denominator for floating point
    real(MK),parameter :: MINDENOM=1e-15_MK
    real(MK),parameter :: MIN_EXP_VALUE=-10.0_MK
    !integer:: iunitSeg
    !logical :: bOutputSeg=.false.
contains
    !>
    !> Takes n segments and n control points
    subroutine fUi_SegmentCst()
        implicit none
    end subroutine

    !> Takes n segments one ofter the other and n control points
    subroutine fUi_ContinousLine(CPs,SgmtP,Intensities, SmoothModel,SmoothParam,bComputeGrad,Uind, Grad,nCPs,nSgmt)
        implicit none
        ! Arguments declarations 
        integer, intent(in) :: nCPs
        integer, intent(in) :: nSgmt
        logical, intent(in) :: bComputeGrad !< 
        real(MK), dimension(3,nCPs), intent(in)    :: CPs          !< side effects
        real(MK), dimension(3,nCPs), intent(inout) :: Uind         !< side effects
        real(MK), dimension(9,nCPs), intent(inout) :: Grad         !< side effects
        real(MK), dimension(3,nSgmt+1), intent(in) :: SgmtP        !< 
        real(MK), dimension(nSgmt), intent(in)   :: Intensities  !< 
        integer, intent(in)                        :: SmoothModel  !< 
        real(MK), dimension(nSgmt), intent(in)     :: SmoothParam  !< 
        ! Variable declarations 
        real(MK), dimension(3) :: Uindtmp      !<  
        real(MK), dimension(9) :: Gradtmp      !<  
        real(MK), dimension(3) :: DP1, DP2
        integer :: icp,is
        ! loop on control points 
        do icp=1,nCPs
            ! loop on segments
            do is=1,nSgmt
                Uindtmp=0.0_MK;
                Gradtmp=0.0_MK;
                DP1=CPs(:,icp)-SgmtP(:,is)
                DP2=CPs(:,icp)-SgmtP(:,is+1)
                call fUi_SegmentCst_11(DP1, DP2, Intensities(is), SmoothModel,&
                    SmoothParam(is), bComputeGrad, Uindtmp, Gradtmp)
                Uind(1:3,icp) = Uind(1:3,icp)+Uindtmp(1:3)
                Grad(1:9,icp) = Grad(1:9,icp)+Gradtmp(1:9)
            enddo 
        enddo
    end subroutine


    !> Used to be fUi_VortexLine, takes one segment and one control points
    subroutine fUi_SegmentCst_11_b(DeltaPa, DeltaPb, norm2_r0, Intensity, visc_model , visc_param, bComputeGrad, Ui, Grad)
        use MathConstants 
        implicit none 
        ! Input/output arguments 
        real(MK), dimension(3), intent(in) :: DeltaPa      !< 3 x 1   Pcp-P1 !!!!!!!!!!
        real(MK), dimension(3), intent(in) :: DeltaPb      !< 3 x 1   Pcp-P2
        real(MK), intent(in)               :: Intensity    !< 
        real(MK), intent(in)               :: norm2_r0      !<  norm squared of segment
        integer, intent(in)                :: visc_model   !< 
        real(MK), intent(in)               :: visc_param   !< 
        logical, intent(in)                :: bComputeGrad !< 
        real(MK), dimension(3),intent(out) :: Ui           !<  No Side effects
        real(MK), dimension(9),intent(out) :: Grad         !<  No Side effects
        ! Variables declaration 
        real(MK),dimension(3) :: crossprod      !< 
        !real(MK),dimension(3) :: D              !< 
        real(MK)              :: denom          !< 
        real(MK)              :: denominator    !< 
        real(MK)              :: h2             !< Square of h
!         real(MK)              :: h2_inv         !< Square of h
        real(MK)              :: h              !< Only used by one model
        real(MK)              :: Kv             !< 
        !real(MK)              :: Kv2            !< Gradient kernel only for Lamb-Oseen
!         real(MK)              :: Kv_tmp         !< 
        real(MK)              :: norm_a         !< 
        real(MK)              :: norm_b         !< 
        real(MK)              :: norm2_crossprod !< 
        !real(MK)              :: norm2_r0tmp !< 
        real(MK)              :: xa             !< 
        real(MK)              :: ya             !< 
        real(MK)              :: za             !< 
        real(MK)              :: xb             !< 
        real(MK)              :: yb             !< 
        real(MK)              :: zb             !< 
        real(MK)              :: visc_param2    !< Square of viscous param
        real(MK)              :: exp_value      !< 
        ! 
        xa=DeltaPa(1)
        ya=DeltaPa(2)
        za=DeltaPa(3)
        xb=DeltaPb(1)
        yb=DeltaPb(2)
        zb=DeltaPb(3)
        !--------------------------------------------------------------------------------
        !--- Simple test if on an extremity point
        !--------------------------------------------------------------------------------
        if(abs(xa)+abs(ya)+abs(za)<MINNORMSIMP) then
            Ui(1:3)=0.0_MK
            Grad(1:9)=0.0_MK
        elseif(abs(xb)+abs(yb)+abs(zb)<MINNORMSIMP) then
            Ui(1:3)=0.0_MK
            Grad(1:9)=0.0_MK
        else
            norm_b = sqrt(xb*xb + yb*yb + zb*zb)
            norm_a = sqrt(xa**2 + ya**2 + za**2)
            !denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
            denom       =               (norm_a*norm_b + xa*xb+ya*yb+za*zb)
            denominator = norm_a*norm_b*denom
            crossprod(1) = ya*zb-za*yb
            crossprod(2) = za*xb-xa*zb
            crossprod(3) = xa*yb-ya*xb
            norm2_crossprod=crossprod(1)**2+ crossprod(2)**2+ crossprod(3)**2
            !print "(A,6F21.16)","xa xb:",xa,ya,za,xb,yb,zb
            ! check for singularity  
            ! TODO this check is viscous model dependent(e.g. crossprod)... have to go towards different funcitons soon..
            ! to avoid singularity, save a bit of computation and use limit values
            if (denominator<PRECISION_UI .or. norm2_crossprod<PRECISION_UI) then
                !--------------------------------------------------------------------------------
                !--- Exactly on the Singularity, velocity is zero, but gradient may be constant (is is wanted?)
                !--------------------------------------------------------------------------------
                !                 if(denominator<MINNORM) then
                !                     print*,'because of denom'
                !                 endif
                !                 if(norm2_crossprod<MINNORM2) then
                !                     print*,'because of crossprod'
                !                 endif
                Ui(1:3)=0.0_MK
                Grad(1:9)=0.0_MK
                if( visc_model==0) then 
                    return
                end if
                if(bComputeGrad) then
                    ! denom=(norm_a*norm_b + xa*xb+ya*yb+za*zb)+MINDENOM
                    ! if(norm_a <MINNORM) norm_a=MINNORM
                    ! if(norm_b <MINNORM) norm_b=MINNORM
                    ! if(denom <MINNORM) denom=MINNORM
                    ! ! For now only one limit for all models...
                    ! h2 = MINNORM2
                    ! Kv = 1.0_MK-exp(-1.25643_MK*h2/visc_param2)
                    ! Grad(2) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(zb-za) ! 12
                    ! Grad(3) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(ya-yb) ! 13
                    ! Grad(4) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(za-zb) ! 21
                    ! Grad(6) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(xb-xa) ! 23
                    ! Grad(7) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(yb-ya) ! 31
                    ! Grad(8) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(xa-xb) ! 32
                    ! print*,Grad(4)
                endif
            else
                !norm2_r0tmp = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb)
                ! if(abs(norm2_r0-norm2_r0tmp)>1e-3) then
                !     print*,'Norm mismatch',norm2_r0tmp,norm2_r0
                !     STOP
                ! endif
                denominator=denominator+MINDENOM
                ! TODO far distance
                !--------------------------------------------------------------------------------
                !--- Normal Procedure 
                !--------------------------------------------------------------------------------
                ! viscous model  
                select case ((visc_model)) !
                case ( 0 ) !! No vortex core model
                    Kv=1.0_MK
                case ( 1 ) !!Rankine - t<=>rc
                    ! orthogonal distance r1xr2/r0 
                    h2 = norm2_crossprod/ norm2_r0 
                    visc_param2=visc_param**2
                    if (h2<visc_param2) then 
                        Kv=h2/visc_param2
                    else
                        Kv=1.0_MK 
                    end if 
                case ( 2 ) !!Lamb-Oseen - vsic_param<=>rc 
                    ! BUG WITH OPENMP ??
                    ! orthogonal distance r1xr2/r0 
                    h2 = norm2_crossprod/ norm2_r0 
                    visc_param2=visc_param**2
                    exp_value = -1.25643_MK*(h2)/(visc_param2)
                    if(exp_value<MIN_EXP_VALUE) then
                        Kv = 1.0_MK
                    else
                        Kv = 1.0_MK-exp(exp_value)
                    endif
                case ( 3 ) !!Vatistas n=2 - visc_param<=>rc 
                    ! orthogonal distance r1xr2/r0 
                    h2 = norm2_crossprod/ norm2_r0 
                    visc_param2=visc_param**2
                    ! h = (norm_a+norm_b)/2; 
                    Kv = h2/sqrt(visc_param2**2+h2**2)
                case ( 4 ) !!Cut-off radius as function of length: delta^2*norm(r0)^2  
                    Kv=1.0_MK
                    denominator=denominator+visc_param*visc_param*norm2_r0
                case ( 5 ) !!Cut-off radius dimension fixed: (delta l_0)^2  
                    Kv=1.0_MK
                    denominator=denominator+visc_param*visc_param
                case ( 33 ) !!Vatistas n=2 - visc_paramt<=>rc  
                    ! See Hoydonck 2012. and Script MainRingSensitivityN
                    ! Not matture enough
                    h = (norm_a+norm_b)/2.0_MK
                    visc_param2=visc_param**2
                    if(h<2*sqrt(norm2_r0)) then
                        ! use Vatistas n=2, normal one
                        h2 = norm2_crossprod/norm2_r0
                    end if
                    Kv = (h2)/sqrt(visc_param2**2+h2**2)
                case DEFAULT
                    Kv=1.0_MK
                end select 

                Kv=Intensity*fourpi_inv*Kv*(norm_a+norm_b)/denominator
                Ui(1:3) = Kv*crossprod(1:3)
                !print "(A,3F21.16)","ui   :",Ui(1:3)

                if (bComputeGrad) then 
                    Grad(1:9)=0.0_MK
                    !denom=(norm_a*norm_b + xa*xb+ya*yb+za*zb)+MINDENOM
!                     denom=denom+MINDENOM
!                     D=-(DeltaPa/norm_a**3+DeltaPb/norm_b**3)&
!                         -(1./norm_a+1./norm_b)*1/denom*(DeltaPa/norm_a*norm_b+ DeltaPb/norm_b*norm_a+ DeltaPa+DeltaPb)
! 
!                     Grad(1) = crossprod(1)*D(1)                              ! 11
!                     Grad(2) = crossprod(1)*D(2) +(1./norm_a+1./norm_b)*(zb-za) ! 12
!                     Grad(3) = crossprod(1)*D(3) +(1./norm_a+1./norm_b)*(ya-yb) ! 13
! 
!                     Grad(4) = crossprod(2)*D(1) +(1./norm_a+1./norm_b)*(za-zb) ! 21
!                     Grad(5) = crossprod(2)*D(2)                              ! 22
!                     Grad(6) = crossprod(2)*D(3) +(1./norm_a+1./norm_b)*(xb-xa) ! 23
! 
!                     Grad(7) = crossprod(3)*D(1) +(1./norm_a+1./norm_b)*(yb-ya) ! 31
!                     Grad(8) = crossprod(3)*D(2) +(1./norm_a+1./norm_b)*(xa-xb) ! 32
!                     Grad(9) = crossprod(3)*D(3)                              ! 33
! 
!                     Grad=Intensity*fourpi_inv*1./denom*Grad;
!                     if(visc_model>0) then
!                         ! h2 and visc_param2 are the same for everybody
!                         ! We use Lamb-Oseen model (it's the easiest)
!                         exp_value = -1.25643_MK*(h2)/(visc_param2)
!                         if(exp_value<MIN_EXP_VALUE) then
!                             ! Far away, we then use inviscid Gradient from above.
!                         else
!                             ! The Lamb-Oseen kernel (may be different than above Kv)
!                             Kv2 = 1.0_MK-exp(exp_value)
!                             ! Computing the Grad of Kv2
!                             D=0._MK
!                             ! Grad(raxrb) ra x rb 
!                             ! !!!!! Watcth out, need the segment extremity coordinates, but we don't have it
!                             ! It's the opposite formula hence, since xa-xb=x1-x2
!                             D(1)=(zb-za)*crossprod(2)+(ya-yb)*crossprod(3);
!                             D(2)=(za-zb)*crossprod(1)+(xb-xa)*crossprod(3);
!                             D(3)=(yb-ya)*crossprod(1)+(xa-xb)*crossprod(2);
!                             D=-2._MK*1.25643_MK/(norm2_r0 *visc_param2)*exp(exp_value)*D;
!                             !                         print*,'Kv',Kv
!                             !                         print*,'Grad',Grad
! 
!                             ! Grad_viscous=Kv2*Grad + u_invisc * D
!                             Grad=Kv2*Grad ! Intensity already in Grad and Ui
                            ! Since Ui contains already Kv, we have to divide by it. Alternative: store Ui inviscid
!                            print*,'Kv',Kv,visc_param,exp_value,Kv2! this print is here to remind that c code has no grad
                           ! also a singulartiy can occur with the grid.
!                             Grad(1) = Grad(1) + Ui(1)*D(1)/Kv
!                             Grad(2) = Grad(2) + Ui(1)*D(2)/Kv
!                             Grad(3) = Grad(3) + Ui(1)*D(3)/Kv
! 
!                             Grad(4) = Grad(4) + Ui(2)*D(1)/Kv
!                             !                         print*, 'Term3', Ui(2)*D(1)/Kv
!                             Grad(5) = Grad(5) + Ui(2)*D(2)/Kv
!                             Grad(6) = Grad(6) + Ui(2)*D(3)/Kv
! 
!                             Grad(7) = Grad(7) + Ui(3)*D(1)/Kv
!                             Grad(8) = Grad(8) + Ui(3)*D(2)/Kv
!                             Grad(9) = Grad(9) + Ui(3)*D(3)/Kv
!                         endif
!                     end if  ! Visc model >0
                else
                    Grad(1:9)=0.0_MK
                end if  ! bComputeGrad
            end if ! denominator size or distances too small
        end if ! 
    end subroutine fUi_SegmentCst_11_b

    !> Used to be fUi_VortexLine, takes one segment and one control points
    subroutine fUi_SegmentCst_11(DeltaPa, DeltaPb, Intensity, visc_model , visc_param, bComputeGrad, Ui, Grad)
        use MathConstants 
        implicit none 
        ! Input/output arguments 
        real(MK), dimension(3), intent(in) :: DeltaPa      !< 3 x 1   Pcp-P1 !!!!!!!!!!
        real(MK), dimension(3), intent(in) :: DeltaPb      !< 3 x 1   Pcp-P2
        real(MK), intent(in)               :: Intensity    !< 
        integer, intent(in)                :: visc_model   !< 
        real(MK), intent(in)               :: visc_param   !< 
        logical, intent(in)                :: bComputeGrad !< 
        real(MK), dimension(3),intent(out) :: Ui           !<  No Side effects
        real(MK), dimension(9),intent(out) :: Grad         !<  No Side effects
        ! Variables declaration 
        real(MK),dimension(3) :: crossprod      !< 
        real(MK),dimension(3) :: D              !< 
        real(MK)              :: denom          !< 
        real(MK)              :: denominator    !< 
        real(MK)              :: h2             !< Square of h
        real(MK)              :: h              !< Only used by one model
        real(MK)              :: Kv             !< 
        real(MK)              :: norm_a         !< 
        real(MK)              :: norm_b         !< 
        real(MK)              :: norm2_r0        !< 
        real(MK)              :: norm2_crossprod !< 
        real(MK)              :: xa             !< 
        real(MK)              :: ya             !< 
        real(MK)              :: za             !< 
        real(MK)              :: xb             !< 
        real(MK)              :: yb             !< 
        real(MK)              :: zb             !< 
        real(MK)              :: visc_param2    !< Square of viscous param
        real(MK)              :: exp_value      !< 
        ! 
        xa=DeltaPa(1)
        ya=DeltaPa(2)
        za=DeltaPa(3)
        xb=DeltaPb(1)
        yb=DeltaPb(2)
        zb=DeltaPb(3)
        !--------------------------------------------------------------------------------
        !--- Simple test if on an extremity point
        !--------------------------------------------------------------------------------
        if(abs(xa)+abs(ya)+abs(za)<MINNORMSIMP) then
            Ui(1:3)=0.0_MK
            Grad(1:9)=0.0_MK
            !print*,'Exit1'
        elseif(abs(xb)+abs(yb)+abs(zb)<MINNORMSIMP) then
            Ui(1:3)=0.0_MK
            Grad(1:9)=0.0_MK
            !print*,'Exit2'
        else
            norm_b = sqrt(xb*xb + yb*yb + zb*zb)
            norm_a = sqrt(xa**2 + ya**2 + za**2)
            norm2_r0 = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb)
            !denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
            denom       =               (norm_a*norm_b + xa*xb+ya*yb+za*zb)
            denominator = norm_a*norm_b*denom
            crossprod(1) = ya*zb-za*yb
            crossprod(2) = za*xb-xa*zb
            crossprod(3) = xa*yb-ya*xb
            norm2_crossprod=crossprod(1)**2+ crossprod(2)**2+ crossprod(3)**2
            !if(abs(norm2_r0-norm2_seg)>1e-20_MK) then
            !    print*,'Norm mismatch'
            !endif
            ! check for singularity  
            !if(bOutputSeg) then
            !    write(iunitSeg,'(a,3EN14.2E2)') 'a', xa,ya,za
            !    write(iunitSeg,'(a,3EN14.2E2)') 'b', xb,yb,zb
            !    write(iunitSeg,'(a,3EN14.2E2)') 'sum', norm_a+norm_b, sqrt(norm2_r0),norm_a+norm_b-sqrt(norm2_r0)
            !    write(iunitSeg,'(3EN14.2E2)') norm2_r0 , norm_a, norm_b
            !    write(iunitSeg,'(a, 2EN14.2E2)') 'den', denominator, MINNORM*norm2_r0 
            !    write(iunitSeg,'(a, 3EN14.2E2)') 'cor', norm2_crossprod, MINNORM2*norm2_r0, PRECISION_UI
            !endif
            ! TODO this check is viscous model dependent(e.g. crossprod)... have to go towards different funcitons soon..
            ! to avoid singularity, save a bit of computation and use limit values
            if (norm2_r0<PRECISION_UI .or. denominator<PRECISION_UI .or. norm2_crossprod<PRECISION_UI) then
                !--------------------------------------------------------------------------------
                !--- Exactly on the Singularity, velocity is zero, but gradient may be constant (is is wanted?)
                !--------------------------------------------------------------------------------
                if(norm2_r0<PRECISION_UI) then
                    !print*,'Exit3, because of segment too small'
                else
                    if(denominator<PRECISION_UI) then
                        !print*,'Exit3, because of denom'
                    endif
                    if(norm2_crossprod<PRECISION_UI) then
                        !print*,'Exit3, because of crossprod'
                    endif
                endif
                !if(bOutputSeg) then
                !    write(iunitSeg,'(a)') 'exitc'
                !endif
                Ui(1:3)=0.0_MK
                Grad(1:9)=0.0_MK
                if( visc_model==0) then 
                    return
                end if
                if(bComputeGrad) then
                    ! denom=(norm_a*norm_b + xa*xb+ya*yb+za*zb)+MINDENOM
                    ! if(norm_a <MINNORM) norm_a=MINNORM
                    ! if(norm_b <MINNORM) norm_b=MINNORM
                    ! if(denom <MINNORM) denom=MINNORM
                    ! ! For now only one limit for all models...
                    ! h2 = MINNORM2
                    ! Kv = 1.0_MK-exp(-1.25643_MK*h2/visc_param2)
                    ! Grad(2) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(zb-za) ! 12
                    ! Grad(3) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(ya-yb) ! 13
                    ! Grad(4) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(za-zb) ! 21
                    ! Grad(6) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(xb-xa) ! 23
                    ! Grad(7) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(yb-ya) ! 31
                    ! Grad(8) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(xa-xb) ! 32
                    ! print*,Grad(4)
                endif
            else
                !denominator=denominator+MINDENOM
                denominator=denominator
                ! TODO far distance
                !--------------------------------------------------------------------------------
                !--- Normal Procedure 
                !--------------------------------------------------------------------------------
                ! viscous model  
                select case ((visc_model)) !
                case ( 0 ) !! No vortex core model
                    Kv=1.0_MK
                case ( 1 ) !!Rankine - t<=>rc
                    ! orthogonal distance r1xr2/r0 
                    h2 = norm2_crossprod/ norm2_r0 
                    visc_param2=visc_param**2
                    if (h2<visc_param2) then 
                        Kv=h2/visc_param2
                    else
                        Kv=1.0_MK 
                    end if 
                case ( 2 ) !!Lamb-Oseen - vsic_param<=>rc 
                    ! BUG WITH OPENMP ??
                    ! orthogonal distance r1xr2/r0 
                    h2 = norm2_crossprod/ norm2_r0 
                    visc_param2=visc_param**2
                    exp_value = -1.25643_MK*(h2)/(visc_param2)
                    if(exp_value<MIN_EXP_VALUE) then
                        Kv = 1.0_MK
                    else
                        Kv = 1.0_MK-exp(exp_value)
                    endif
                case ( 3 ) !!Vatistas n=2 - visc_param<=>rc 
                    ! orthogonal distance r1xr2/r0 
                    h2 = norm2_crossprod/ norm2_r0 
                    visc_param2=visc_param**2
                    ! h = (norm_a+norm_b)/2; 
                    Kv = h2/sqrt(visc_param2**2+h2**2)
                case ( 4 ) !!Cut-off radius no dimensions - visc_param<=>delta**2 
                    h2 = norm2_crossprod/ norm2_r0 
                    Kv=1.0_MK
                    ! delta*norm(r0)^2  
                    denominator=denominator+visc_param*norm2_r0
                    visc_param2=visc_param**2 ! TODO
                case ( 5 ) !!Cut-off radius dimension - visc_param<=>(delta l)**2 
                    h2 = norm2_crossprod/ norm2_r0 
                    Kv=1.0_MK
                    ! (delta l_0)^2  
                    denominator=denominator+visc_param
                    visc_param2=visc_param**2 ! TODO
                case ( 33 ) !!Vatistas n=2 - visc_paramt<=>rc  
                    ! See Hoydonck 2012. and Script MainRingSensitivityN
                    ! Not matture enough
                    h = (norm_a+norm_b)/2.0_MK
                    visc_param2=visc_param**2
                    if(h<2*sqrt(norm2_r0)) then
                        ! use Vatistas n=2, normal one
                        h2 = norm2_crossprod/norm2_r0
                    else
                        h2 = 1._MK ! TODO
                    endif
                    Kv = (h2)/sqrt(visc_param2**2+h2**2)
                case DEFAULT
                    Kv=1.0_MK
                end select 

                Kv=Intensity*fourpi_inv*Kv*(norm_a+norm_b)/denominator
                Ui(1:3) = Kv*crossprod(1:3)
                !if(bOutputSeg) then
                !    write(iunitSeg,'(A,E14.2E2)') 'Kv', Kv
                !    write(iunitSeg,'(A,3E14.2E2)') 'crossprod', crossprod
                !    write(iunitSeg,'(A,E14.2E2)') 'Intensity', Intensity
                !    write(iunitSeg,'(A,3E14.2E2)') 'Ui', Ui
                !endif

                if (bComputeGrad) then 
                    denom=(norm_a*norm_b + xa*xb+ya*yb+za*zb)+MINDENOM
                    D=-(DeltaPa/norm_a**3+DeltaPb/norm_b**3)&
                        -(1/norm_a+1/norm_b)*1/denom*(DeltaPa/norm_a*norm_b+ DeltaPb/norm_b*norm_a+ DeltaPa+DeltaPb)

                    Grad(1) = crossprod(1)*D(1)                              ! 11
                    Grad(2) = crossprod(1)*D(2) +(1/norm_a+1/norm_b)*(zb-za) ! 12
                    Grad(3) = crossprod(1)*D(3) +(1/norm_a+1/norm_b)*(ya-yb) ! 13

                    Grad(4) = crossprod(2)*D(1) +(1/norm_a+1/norm_b)*(za-zb) ! 21
                    Grad(5) = crossprod(2)*D(2)                              ! 22
                    Grad(6) = crossprod(2)*D(3) +(1/norm_a+1/norm_b)*(xb-xa) ! 23

                    Grad(7) = crossprod(3)*D(1) +(1/norm_a+1/norm_b)*(yb-ya) ! 31
                    Grad(8) = crossprod(3)*D(2) +(1/norm_a+1/norm_b)*(xa-xb) ! 32
                    Grad(9) = crossprod(3)*D(3)                              ! 33

                    Grad=Intensity*fourpi_inv*1/denom*Grad;
                    if(visc_model>0) then
                        ! h2 and visc_param2 are the same for everybody
                        ! We use Lamb Oseen model !!! TODO, expensive if far away!!!
                        exp_value = -1.25643_MK*(h2)/(visc_param2)
                        if(exp_value<MIN_EXP_VALUE) then
                            ! Far away, not need to compute smooth gradient
                            Kv = 1.0_MK
                        else
                            Kv = 1.0_MK-exp(exp_value)
                        !                     end select
                            ! Computing the Grad of Kv
                            D=0._MK
                            ! Grad(raxrb) ra x rb 
                            ! !!!!! Watcth out, need the segment extremity coordinates, but we don't have it
                            ! It's the opposite formula hence, since xa-xb=x1-x2
                            D(1)=(zb-za)*crossprod(2)+(ya-yb)*crossprod(3);
                            D(2)=(za-zb)*crossprod(1)+(xb-xa)*crossprod(3);
                            D(3)=(yb-ya)*crossprod(1)+(xa-xb)*crossprod(2);
                            D=-2*1.25643/(norm2_r0 *visc_param2)*exp(exp_value)*D;
                            !                         print*,'Kv',Kv
                            !                         print*,'Grad',Grad

                            Grad=Kv*Grad ! Intensity already in Grad and Ui
                            ! Since Ui contains already Kv, we have to divide by it..
                            Grad(1) = Grad(1) + Ui(1)*D(1)/Kv
                            Grad(2) = Grad(2) + Ui(1)*D(2)/Kv
                            Grad(3) = Grad(3) + Ui(1)*D(3)/Kv

                            Grad(4) = Grad(4) + Ui(2)*D(1)/Kv
                            !                         print*, 'Term3', Ui(2)*D(1)/Kv
                            Grad(5) = Grad(5) + Ui(2)*D(2)/Kv
                            Grad(6) = Grad(6) + Ui(2)*D(3)/Kv

                            Grad(7) = Grad(7) + Ui(3)*D(1)/Kv
                            Grad(8) = Grad(8) + Ui(3)*D(2)/Kv
                            Grad(9) = Grad(9) + Ui(3)*D(3)/Kv
                        endif
                    end if  ! Visc model >0
                else
                    Grad(1:9)=0.0_MK
                end if  ! bComputeGrad
            end if ! denominator size or distances too small
        end if ! 
        ! printf("!4.3f !4.3f !4.3f !4.3f !4.3fNewLine",Uout(1),Uout(2),Uout(3),Kv,denominator);  
    end subroutine fUi_SegmentCst_11


end module UISegments




<|start_filename|>welib/vortilib/elements/fortran/UIPointSource.f90<|end_filename|>
!> 
module UIPointSource
    use PrecisionMod, only: MK
    implicit none

    real(MK),parameter :: MINNORM=1e-4
contains

    subroutine fUi_PointSource(CPs,SrcP,Intensities, SmoothModel,SmoothParam,bComputeGrad,Uind, Grad,nCPs,nSrc)
        ! Arguments declarations 
        integer, intent(in) :: nCPs
        integer, intent(in) :: nSrc
        logical, intent(in) :: bComputeGrad !< 
        real(MK), dimension(3,nCPs), intent(in)    :: CPs         !< side effects
        real(MK), dimension(3,nCPs), intent(inout) :: Uind        !< side effects
        real(MK), dimension(9,nCPs), intent(inout) :: Grad        !< side effects
        real(MK), dimension(3,nSrc+1), intent(in)  :: SrcP        !< 
        real(MK), dimension(nSrc), intent(in)      :: Intensities !< 
        integer, intent(in)                        :: SmoothModel !< 
        real(MK), dimension(nSrc), intent(in)      :: SmoothParam !< 
        ! Variable declarations 
        integer :: icp,is
        real(MK), dimension(3) :: Uindtmp !< 
        real(MK), dimension(9) :: Gradtmp !< 
        real(MK), dimension(3) :: DP      !< 
        ! loop on control points 
        do icp=1,nCPs
            ! loop on Sources
            do is=1,nSrc
                Uindtmp=0.0_MK;
                Gradtmp=0.0_MK;
                DP=CPs(1:3,icp)-SrcP(1:3,is)
                call fUi_PointSource_11(DP,Intensities(is),SmoothModel,SmoothParam(is),bComputeGrad,&
                    Uindtmp,Gradtmp)
                Uind(1:3,icp) = Uind(1:3,icp)+Uindtmp(1:3)
                Grad(1:9,icp) = Grad(1:9,icp)+Gradtmp(1:9)
            enddo 
        enddo
    end subroutine

    subroutine fUi_PointSource_11(DeltaP,Intensity, SmoothModel,SmoothParam,bComputeGrad,Uind, Grad)
        use MathConstants ,only: fourpi_inv
        ! Input/output arguments 
        real(MK), dimension(3), intent(in) :: DeltaP       !< 3 x 1   Pcp-Ps
        real(MK), intent(in)               :: Intensity    !< 
        integer, intent(in)                :: SmoothModel  !< 
        real(MK), intent(in)               :: SmoothParam  !< 
        logical, intent(in)                :: bComputeGrad !< 
        real(MK), dimension(3),intent(out) :: Uind         !< No Side effects
        real(MK), dimension(9),intent(out) :: Grad         !< No Side effects
        ! Variables declaration 
        real(MK) :: r      !< 
        real(MK) :: r3_inv !< 
        real(MK) :: r2_inv !< 
        real(MK) :: E      !< 
        real(MK) :: Ebar   !< 
        real(MK) :: factor   !< 
        real(MK) :: Kv   !< Mollifier

        r=sqrt(DeltaP(1)**2+ DeltaP(2)**2+ DeltaP(3)**2)

        Uind(1:3)=0.0_MK
        Grad(1:9)=0.0_MK
        ! TODO: escape should depend on physics
        if (r<MINNORM) then! to avoid singularity, save a bit of computation and use limit values
            return
            !--------------------------------------------------------------------------------
            !--- Exactly on the Singularity 
            !--------------------------------------------------------------------------------
            !if( SmoothModel==0) then 
            !    return
            !end if
            ! NOTE: We use to return the value obtained analytically for the gradient on the singularity
            ! Now, we just return 0
            !
            !if(bComputeGrad) then
            !    !if( SmoothModel==2) then 
            !    !    r3_inv=1._MK/SmoothParam**3 ! Epsilon^3
            !    !else
            !    !    r3_inv=1._MK/SmoothParam**3 ! Epsilon ^3
            !    !endif
            !    ! For model 1 and 2 is the same limit
            !    !Grad(1) = fourpi_inv*Intensity*r3_inv ! 11
            !    !Grad(5) = fourpi_inv*Intensity*r3_inv ! 22
            !    !Grad(9) = fourpi_inv*Intensity*r3_inv ! 33
            !endif
            ! TODO far distance
        else
            !--------------------------------------------------------------------------------
            !--- Normal Procedure 
            !--------------------------------------------------------------------------------
            ! Velocity 
            select case(SmoothModel)
            case(0)
                r3_inv=1._MK/(r**3);
                Kv=1.0_MK
            case (1) !
                !Exponential mollifier SmoothParam=Epsilon^3
                r3_inv=1._MK/(r**3)
                E=exp(-r**3/SmoothParam**3)
                Ebar=1._MK-E
                Kv=Ebar
            case (2) !
                ! Finite support with SmoothParam=Epsilon^6  !!!
                r3_inv= 1._MK/sqrt(SmoothParam**6+r**6)
                Kv= 1.0_MK
            case default
                print*,'Wrong smooth model for point source'
            end select

            Uind(1:3)=Kv*Intensity*fourpi_inv*(DeltaP)*r3_inv 

            ! Gradient 
            if (bComputeGrad) then 
                select case ( SmoothModel) !
                case (0) !
                    r2_inv=1._MK/r**2
                    Grad(1) = fourpi_inv*r3_inv*Intensity*(1 -3*DeltaP(1)*DeltaP(1)*r2_inv) ! 11
                    Grad(2) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(2)*DeltaP(1)*r2_inv) ! 12
                    Grad(3) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(3)*DeltaP(1)*r2_inv) ! 13

                    Grad(4) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(1)*DeltaP(2)*r2_inv) ! 21
                    Grad(5) = fourpi_inv*r3_inv*Intensity*(1 -3*DeltaP(2)*DeltaP(2)*r2_inv) ! 22
                    Grad(6) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(3)*DeltaP(2)*r2_inv) ! 23

                    Grad(7) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(1)*DeltaP(3)*r2_inv) ! 31
                    Grad(8) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(2)*DeltaP(3)*r2_inv) ! 32
                    Grad(9) = fourpi_inv*r3_inv*Intensity*(1 -3*DeltaP(3)*DeltaP(3)*r2_inv) ! 33

                case (1) !
                    ! Exponential mollifier 
                    r2_inv=1._MK/r**2
                    ! Almost the same as above 
                    ! Remember I chose SmoothParam=epsilon**3, so don't cube it twice!!
                    factor=-3*Ebar*r2_inv+3*r*E/SmoothParam**3
                    Grad(1) = fourpi_inv*r3_inv*Intensity*( Ebar  + factor*DeltaP(1)*DeltaP(1))!11
                    Grad(2) = fourpi_inv*r3_inv*Intensity*(         factor*DeltaP(2)*DeltaP(1))!12
                    Grad(3) = fourpi_inv*r3_inv*Intensity*(         factor*DeltaP(3)*DeltaP(1))!13

                    Grad(4) = fourpi_inv*r3_inv*Intensity*(         factor*DeltaP(1)*DeltaP(2))!21
                    Grad(5) = fourpi_inv*r3_inv*Intensity*( Ebar  + factor*DeltaP(2)*DeltaP(2))!22
                    Grad(6) = fourpi_inv*r3_inv*Intensity*(         factor*DeltaP(3)*DeltaP(2))!23

                    Grad(7) = fourpi_inv*r3_inv*Intensity*(         factor*DeltaP(1)*DeltaP(3))!31
                    Grad(8) = fourpi_inv*r3_inv*Intensity*(         factor*DeltaP(2)*DeltaP(3))!32
                    Grad(9) = fourpi_inv*r3_inv*Intensity*( Ebar  + factor*DeltaP(3)*DeltaP(3))!33

                case (2) !
                    ! Finite support with Epsilon 6 
                    r2_inv= r**4/(SmoothParam**6+r**6) ! Fake r2_inv
                    ! EXACTLY the same as case 0 due to fake variables r3 and r5 
                    Grad(1) = fourpi_inv*r3_inv*Intensity*(1 -3*DeltaP(1)*DeltaP(1)*r2_inv) ! 11
                    Grad(2) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(2)*DeltaP(1)*r2_inv) ! 12
                    Grad(3) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(3)*DeltaP(1)*r2_inv) ! 13

                    Grad(4) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(1)*DeltaP(2)*r2_inv) ! 21
                    Grad(5) = fourpi_inv*r3_inv*Intensity*(1 -3*DeltaP(2)*DeltaP(2)*r2_inv) ! 22
                    Grad(6) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(3)*DeltaP(2)*r2_inv) ! 23

                    Grad(7) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(1)*DeltaP(3)*r2_inv) ! 31
                    Grad(8) = fourpi_inv*r3_inv*Intensity*(  -3*DeltaP(2)*DeltaP(3)*r2_inv) ! 32
                    Grad(9) = fourpi_inv*r3_inv*Intensity*(1 -3*DeltaP(3)*DeltaP(3)*r2_inv) ! 33
                end select  ! Smooth Model
            end if  ! bComputeGrad
        end if  ! Norm Big enough
    end subroutine
end module


<|start_filename|>welib/tools/external/Makefile<|end_filename|>
all: test

test:
	python lic_demo.py


install:
	python setup.py build_ext --inplace



<|start_filename|>welib/vortilib/elements/fortran/UIVortexCylinders.f90<|end_filename|>
!> 
module UIVortexCylinders
    use PrecisionMod, only: MK
    implicit none

    real(MK),parameter :: MINNORM=1e-4
contains

    !> Semi Infinite vortex cylinder with tangential vorticiy
    subroutine fUi_VortexCylinderTangSemiInf_11(CP,PCenter,R,gamma_t,N, SmoothModel,SmoothParam,bComputeGrad,Uind,Grad)
        ! Arguments declarations 
        real(MK), dimension(3), intent(in)    :: CP           !< 
        real(MK), dimension(3), intent(in)    :: PCenter      !< PCenter : Position of ring center
        real(MK), intent(in)                  :: R            !< Cylinder radius
        real(MK), intent(in)                  :: gamma_t      !< Tangential vorticity Intensity
        real(MK), dimension(3), intent(in)    :: N            !< Vector norm is gamma_t, direction is normal to the
        integer, intent(in)                   :: SmoothModel  !< 
        real(MK), intent(in)                  :: SmoothParam  !< 
        logical, intent(in)                   :: bComputeGrad !< 
        real(MK), dimension(3), intent(inout) :: Uind         !< side effects
        real(MK), dimension(9), intent(inout) :: Grad         !< side effects
        ! Variable declarations 
        real(MK), dimension(3) :: ex  !<  
        real(MK), dimension(3) :: ey  !< 
        real(MK), dimension(3) :: ez  !<  
        real(MK), dimension(3,3) :: T_r2g
        real(MK), dimension(3,3) :: T_g2r
        integer :: imin
        real(MK), dimension(3) :: CP_loc  !< 
        real(MK) :: R_cp
        real(MK) :: Z_cp
        real(MK) :: phi
        real(MK), dimension(3) :: Uindtmp         !< 
        real(MK), dimension(9) :: Gradtmp         !<
        real(MK) :: urad
        real(MK) :: upsi
        !         ! Note: theories return only vr (index 1) and vz (index 3), there is some cartesian handling though, but does it work ? Is it consistent with the rest of the codes. Does it make use of the inputed vpsi to reduce the error in the angle estimate. 
        ! 
        !! Setting up an orthonornal cartesian system normal to the ring
        ez(1:3)=N(1:3)/norm2(N) ! just in case N is not unitary
        imin=minloc(abs(ez),1)
        ey=ez
        ey(imin)=0
        ey=ey([ 2,3,1 ] )
        ey(imin)=-ey(imin)
        ey=ey/norm2(ey)
        ex=cross(ey,ez)
!         print*,ex
!         print*,ey
!         print*,ez
        !! Now setting up rotation matrix 
        T_r2g(1:3,1)=ex(1:3)
        T_r2g(1:3,2)=ey(1:3)
        T_r2g(1:3,3)=ez(1:3)
        T_g2r=transpose(T_r2g)

        !! Coordinate transformation to a cartesian frame normal to the ring
        CP_loc(1)=T_g2r(1,1)*(CP(1)-PCenter(1)) + T_g2r(1,2)*(CP(2)-PCenter(2)) + T_g2r(1,3)*(CP(3)-PCenter(3))
        CP_loc(2)=T_g2r(2,1)*(CP(1)-PCenter(1)) + T_g2r(2,2)*(CP(2)-PCenter(2)) + T_g2r(2,3)*(CP(3)-PCenter(3))
        CP_loc(3)=T_g2r(3,1)*(CP(1)-PCenter(1)) + T_g2r(3,2)*(CP(2)-PCenter(2)) + T_g2r(3,3)*(CP(3)-PCenter(3))

        !! Coordinate transform to a canonical polar coordinates
        Z_cp=CP_loc(3) ! Since we are already in the ring coordinate system 
        R_cp=sqrt(CP_loc(1)**2+CP_loc(2)**2)
        !print*,'CP_loc',CP_loc

        !! Computing velocity in polar coordinates
        Uindtmp=0.0_MK
        Gradtmp=0.0_MK
        !print*,'gamma_t',gamma_t
        ! We use the Smooth Param of the cylinder no matter the SmoothModel
        if(SmoothModel==0) then
            call fUi_CylinderTangSemiInf11_loc(R_cp, Z_cp, gamma_t, R, SmoothParam, Uindtmp)
            !print*,'Uind pol',Uindtmp
        else
            call fUi_CylinderTangSemiInf11_loc(R_cp, Z_cp, gamma_t, R, SmoothParam, Uindtmp)
        endif

        !! Converting to cartesian 
        phi=atan2(CP_loc(2),CP_loc(1)) 
        !print*,'phi',phi*180/3.141593
        urad=Uindtmp(1)
        upsi=0.0_MK
        Uindtmp(1)=+urad*cos(phi)
        Uindtmp(2)=+urad*sin(phi)
        !print*,'Uind cart',Uindtmp

        !! Going back to global frame 
        Uind(1)=Uind(1) + T_r2g(1,1)*Uindtmp(1) + T_r2g(1,2)*Uindtmp(2) + T_r2g(1,3)*Uindtmp(3)
        Uind(2)=Uind(2) + T_r2g(2,1)*Uindtmp(1) + T_r2g(2,2)*Uindtmp(2) + T_r2g(2,3)*Uindtmp(3)
        Uind(3)=Uind(3) + T_r2g(3,1)*Uindtmp(1) + T_r2g(3,2)*Uindtmp(2) + T_r2g(3,3)*Uindtmp(3)

        !print*,'Uind glob',Uind
        !! TODO gradient 
        if (bComputeGrad) then 
            Grad = 0.0D0  
        else
            Grad = 0.0D0 
        end if 
    contains
        !> Cross product of two vectors ! See crossprod in UTILS...
        function cross(a, b)
            real(MK), dimension(3) :: cross
            real(MK), dimension(3), intent(in) :: a, b
            cross(1) = a(2) * b(3) - a(3) * b(2)
            cross(2) = a(3) * b(1) - a(1) * b(3)
            cross(3) = a(1) * b(2) - a(2) * b(1)
        end function cross
    end subroutine

    !> Induction from a Cylinder at a control point, within the referential of a Cylinder: coordinates (r,theta, z), the ring is at z=0
    subroutine fUi_CylinderTangSemiInf11_loc(Rcp, Zcp, gamma_t, R, SmoothParam, Uind)
        use EllipticIntegrals, only: elliptic_ke, elliptic_pi
        use MathConstants, only:  pi
        ! Arguments declarations 
        real(MK), intent(in) :: Rcp   !< Control point radius
        real(MK), intent(in) :: Zcp   !< Control point z location
        real(MK), intent(in) :: gamma_t !< Cylinder intensity
        real(MK), intent(in) :: R     !< Cylinder radius
        real(MK), intent(in) :: SmoothParam  !< 
        real(MK), dimension(3), intent(inout) :: Uind  !< side effects, in coord (r, theta, z)
        ! Variable declarations 
        real(MK), parameter :: MIN_DIST=1e-7 ! threshold for using axis formula
        real(MK) :: KK  !<  
        real(MK) :: EE  !<  
        real(MK) :: PPI  !<  
        real(MK) :: k  !<  
        real(MK) :: k_2  !<  
        real(MK) :: k0_2  !<  

        if (rcp<MIN_DIST*R) then
            !! Using Axis formula : v_z=-Gamma/(2R) *1 / (1+(z/R)^2)^(3/2) See Naca 1184 - Watsle, or Glauert Contraction of slipstream of an airscrew 1926 
            Uind(1)=0
            Uind(3)=gamma_t/(2_MK)*(1_MK + Zcp/sqrt(Zcp**2+R**2))
        else
            ! Eliptic integrals 
            k_2 = 4_MK*rcp*R/((R+rcp)**2+zcp**2+SmoothParam**2)
            k = sqrt(k_2)
            k0_2 = 4_MK*rcp*R/((R+rcp)**2+SmoothParam**2)
            call elliptic_ke( k_2 , KK,EE) 
            call elliptic_pi(k0_2,90._MK,k_2,PPI)
            ! ur 
            Uind(1)=-gamma_t/(2_MK*pi)*sqrt(R/rcp)*( (2_MK-k_2)/k*KK - 2_MK/k*EE )
            ! ur 
            Uind(3)= gamma_t/2_MK*( 1._MK/2._MK*(1+(R-rcp)/((rcp+R)*sqrt(1-k0_2))) &
                +zcp*k/(2_MK*pi*sqrt(rcp*R))*(KK + (R-rcp)/(R+rcp)*PPI ) )
        endif
    end subroutine fUi_CylinderTangSemiInf11_loc

    subroutine fui_cylinder11_loc_c(Rcp, Zcp, gamma_t, R,SmoothParam, Uind) bind(C,name='fui_cylinder11_loc_c')
        use PrecisionMod, only: C_DOUBLE
        !DEC$ ATTRIBUTES DLLEXPORT :: fui_cylinder11_loc_c
        !GCC$ ATTRIBUTES DLLEXPORT :: fui_cylinder11_loc_c
        ! Arguments declarations 
        real(C_DOUBLE), intent(in) :: Rcp   !< Control point radius
        real(C_DOUBLE), intent(in) :: Zcp   !< Control point z location
        real(C_DOUBLE), intent(in) :: gamma_t !< cylinder intensity
        real(C_DOUBLE), intent(in) :: R     !< cylinder radius
        real(C_DOUBLE), intent(in) :: SmoothParam     !< cylinder radius
        real(C_DOUBLE), dimension(3), intent(inout) :: Uind  !< side effects, in coord (r, theta, z)
        real(MK), dimension(3) :: Uind_MK 

        Uind_MK=0.0_MK
        call fUi_CylinderTangSemiInf11_loc(Rcp, Zcp, gamma_t, R, SmoothParam, Uind_MK)
        Uind=real(Uind_MK,C_DOUBLE)

    end subroutine 

end module


<|start_filename|>Makefile<|end_filename|>
all: test

test:
	python -m unittest discover -v

dep:
	python -m pip install -r requirements.txt

install: dep
	python -m pip install -e .

export:
	python -c 'from welib.tools.repo import *; export_figs_rec("welib/");'



<|start_filename|>welib/vortilib/elements/fortran/Sgmt.f90<|end_filename|>

subroutine fui_sgmt_c(CPs, ioff, iCP_beg, nCPs_loc, nCPs, iComputeGrad, Uind_out, Grad_out, &
        nSgmtTot, nSgmtPTot, SgmtP, Sgmt, iSgmtCurPos, iSgmtStart, &
        SgmtNorm2, SgmtSmoothModel, SgmtSmoothParam, SgmtIntensities)
    use MathConstants 
    use PrecisionMod, only: C_DOUBLE, C_INT, PRECISION_UI
    implicit none
    integer(C_INT), intent(in)                        :: iComputeGrad    !< 
    integer(C_INT), intent(in)                        :: ioff
    integer(C_INT), intent(in)                        :: iCP_beg
    integer(C_INT), intent(in)                        :: nCPs_loc
    integer(C_INT), intent(in)                        :: nCPs
    real(C_DOUBLE), dimension(3,nCPs), intent(in)     :: CPs             !< 
    real(C_DOUBLE), dimension(3,nCPs), intent(inout)  :: Uind_out      !< Induced velocity vector - Side effects!!!
    real(C_DOUBLE), dimension(9,nCPs), intent(inout)  :: Grad_out      !< Gradient vector - Side effects
    integer(C_INT), intent(in)                        :: nSgmtTot
    integer(C_INT), intent(in)                        :: nSgmtPTot
    integer(C_INT), dimension(2,nSgmtTot), intent(in)  :: Sgmt
    real(C_DOUBLE), dimension(3,nSgmtPTot), intent(in) :: SgmtP
    integer(C_INT),intent(in)                         :: iSgmtCurPos
    integer(C_INT),intent(in)                         :: iSgmtStart
    real(C_DOUBLE), dimension(nSgmtTot), intent(in)   :: SgmtNorm2
    integer(C_INT), intent(in)                        :: SgmtSmoothModel
    real(C_DOUBLE), dimension(nSgmtTot), intent(in)   :: SgmtSmoothParam
    real(C_DOUBLE), dimension(nSgmtTot), intent(in)   :: SgmtIntensities
    ! Variables
    integer(C_INT) :: icp,is
    ! --------------------------------------------------------------------------------
    ! --- Variables for inlined C function 
    ! --------------------------------------------------------------------------------
    ! Kindof arguments
    real(C_DOUBLE), dimension(3) :: Uind          !<
    real(C_DOUBLE), dimension(9) :: Grad          !< 
    real(C_DOUBLE), dimension(3) :: P1            !< Point1 of a given segment
    real(C_DOUBLE), dimension(3) :: P2            !< Point2 of a given segment
    real(C_DOUBLE), dimension(3) :: DP1, DP2
    real(C_DOUBLE)               :: norm2_r0
    real(C_DOUBLE)               :: Intensity
    integer(C_INT)               :: visc_model   !< 
    real(C_DOUBLE)               :: visc_param   !< 
    ! Variables declaration 
    real(C_DOUBLE),dimension(3) :: crossprod      !< 
    real(C_DOUBLE)              :: denom          !< 
    real(C_DOUBLE)              :: denominator    !< 
    real(C_DOUBLE)              :: h2             !< Square of h
    real(C_DOUBLE)              :: h              !< Only used by one model
    real(C_DOUBLE)              :: Kv             !< 
    real(C_DOUBLE)              :: norm_a         !< 
    real(C_DOUBLE)              :: norm_b         !< 
    real(C_DOUBLE)              :: norm2_crossprod !< 
    real(C_DOUBLE)              :: xa             !< 
    real(C_DOUBLE)              :: ya             !< 
    real(C_DOUBLE)              :: za             !< 
    real(C_DOUBLE)              :: xb             !< 
    real(C_DOUBLE)              :: yb             !< 
    real(C_DOUBLE)              :: zb             !< 
    real(C_DOUBLE)              :: visc_param2    !< Square of viscous param
    real(C_DOUBLE)              :: exp_value      !< 
    !
    real(MK),parameter :: MIN_EXP_VALUE=-10.0_C_DOUBLE
    real(MK),parameter :: MINDENOM=1e-15_C_DOUBLE
    real(MK),parameter :: MINNORMSIMP=1e-6_C_DOUBLE

    !print*,'Sgmt fortran',(iCP_beg+1),(iCP_beg+1)+nCPS_loc-1,nCPs
    visc_model=SgmtSmoothModel

    !$OMP PARALLEL default(shared)
    !$OMP do private(&
    !$OMP& icp,is,Uind,Grad,P1,P2,DP1,DP2,norm2_r0,Intensity,visc_param,&
    !$OMP& crossprod,denom,denominator,h2,h,Kv,norm_a,norm_b,norm2_crossprod,xa,ya,za,xb,yb,zb,visc_param2,exp_value&
    !$OMP& ) schedule(runtime)
    ! loop on CPs 
    do icp=(iCP_beg+1),(iCP_beg+1)+nCPS_loc-1
        do is=iSgmtStart,iSgmtCurPos ! loop on selected segments 
            ! Kind of arguments for Segment 11
            Uind=0.0_C_DOUBLE
            Grad=0.0_C_DOUBLE 
            P1         = SgmtP(1:3, Sgmt(1,is)) ! Segment extremity points
            P2         = SgmtP(1:3, Sgmt(2,is))
            DP1        = CPs(:,icp)-P1;
            DP2        = CPs(:,icp)-P2
            norm2_r0   = SgmtNorm2(is)
            Intensity  = SgmtIntensities(is)
            visc_param = SgmtSmoothParam(is)
            ! --------------------------------------------------------------------------------
            ! --- inlining of Segment 11 
            ! --------------------------------------------------------------------------------
            ! 
            xa=DP1(1)
            ya=DP1(2)
            za=DP1(3)
            xb=DP2(1)
            yb=DP2(2)
            zb=DP2(3)
            !--------------------------------------------------------------------------------
            !--- Simple test if on an extremity point
            !--------------------------------------------------------------------------------
            if(abs(xa)+abs(ya)+abs(za)<MINNORMSIMP) then
                !Uind(1:3)=0.0_C_DOUBLE
                !Grad(1:9)=0.0_C_DOUBLE
            elseif(abs(xb)+abs(yb)+abs(zb)<MINNORMSIMP) then
                !Uind(1:3)=0.0_C_DOUBLE
                !Grad(1:9)=0.0_C_DOUBLE
            else
                norm_b = sqrt(xb*xb + yb*yb + zb*zb)
                norm_a = sqrt(xa**2 + ya**2 + za**2)
                !denominator = norm_a*norm_b*(norm_a*norm_b + xa*xb+ya*yb+za*zb)
                denom       =               (norm_a*norm_b + xa*xb+ya*yb+za*zb)
                denominator = norm_a*norm_b*denom
                crossprod(1) = ya*zb-za*yb
                crossprod(2) = za*xb-xa*zb
                crossprod(3) = xa*yb-ya*xb
                norm2_crossprod=crossprod(1)**2+ crossprod(2)**2+ crossprod(3)**2
                !print "(A,6F21.16)","xa xb:",xa,ya,za,xb,yb,zb
                ! check for singularity  
                ! TODO this check is viscous model dependent(e.g. crossprod)... have to go towards different funcitons soon..
                ! to avoid singularity, save a bit of computation and use limit values
                if (denominator<PRECISION_UI .or. norm2_crossprod<PRECISION_UI) then
                    !--------------------------------------------------------------------------------
                    !--- Exactly on the Singularity, velocity is zero, but gradient may be constant (is is wanted?)
                    !--------------------------------------------------------------------------------
                    !                 if(denominator<MINNORM) then
                    !                     print*,'because of denom'
                    !                 endif
                    !                 if(norm2_crossprod<MINNORM2) then
                    !                     print*,'because of crossprod'
                    !                 endif
                    if( visc_model==0) then 
                        !return
                    end if
                    if(iComputeGrad==1) then
                        ! denom=(norm_a*norm_b + xa*xb+ya*yb+za*zb)+MINDENOM
                        ! if(norm_a <MINNORM) norm_a=MINNORM
                        ! if(norm_b <MINNORM) norm_b=MINNORM
                        ! if(denom <MINNORM) denom=MINNORM
                        ! ! For now only one limit for all models...
                        ! h2 = MINNORM2
                        ! Kv = 1.0_C_DOUBLE-exp(-1.25643_C_DOUBLE*h2/visc_param2)
                        ! Grad(2) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(zb-za) ! 12
                        ! Grad(3) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(ya-yb) ! 13
                        ! Grad(4) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(za-zb) ! 21
                        ! Grad(6) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(xb-xa) ! 23
                        ! Grad(7) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(yb-ya) ! 31
                        ! Grad(8) = Kv**2*Intensity*fourpi_inv*1/denom*(1/norm_a+1/norm_b)*(xa-xb) ! 32
                        ! print*,Grad(4)
                    endif
                else
                    !norm2_r0tmp = (xa-xb)*(xa-xb) + (ya-yb)*(ya-yb) +(za-zb)*(za-zb)
                    ! if(abs(norm2_r0-norm2_r0tmp)>1e-3) then
                    !     print*,'Norm mismatch',norm2_r0tmp,norm2_r0
                    !     STOP
                    ! endif
                    denominator=denominator+MINDENOM
                    ! TODO far distance
                    !--------------------------------------------------------------------------------
                    !--- Normal Procedure 
                    !--------------------------------------------------------------------------------
                    ! viscous model  
                    h2 = 0.0_MK
                    select case ((visc_model)) !
                    case ( 0 ) !! No vortex core model
                        Kv=1.0_C_DOUBLE
                    case ( 1 ) !!Rankine - t<=>rc
                        ! orthogonal distance r1xr2/r0 
                        h2 = norm2_crossprod/ norm2_r0 
                        visc_param2=visc_param**2
                        if (h2<visc_param2) then 
                            Kv=h2/visc_param2
                        else
                            Kv=1.0_C_DOUBLE 
                        end if 
                    case ( 2 ) !!Lamb-Oseen - vsic_param<=>rc 
                        ! BUG WITH OPENMP ??
                        ! orthogonal distance r1xr2/r0 
                        h2 = norm2_crossprod/ norm2_r0 
                        visc_param2=visc_param**2
                        exp_value = -1.25643_C_DOUBLE*(h2)/(visc_param2)
                        if(exp_value<MIN_EXP_VALUE) then
                            Kv = 1.0_C_DOUBLE
                        else
                            Kv = 1.0_C_DOUBLE-exp(exp_value)
                        endif
                    case ( 3 ) !!Vatistas n=2 - visc_param<=>rc 
                        ! orthogonal distance r1xr2/r0 
                        h2 = norm2_crossprod/ norm2_r0 
                        visc_param2=visc_param**2
                        ! h = (norm_a+norm_b)/2; 
                        Kv = h2/sqrt(visc_param2**2+h2**2)
                    case ( 4 ) !!Cut-off radius as function of length: delta^2*norm(r0)^2  
                        Kv=1.0_C_DOUBLE
                        denominator=denominator+visc_param*visc_param*norm2_r0
                    case ( 5 ) !!Cut-off radius dimension fixed: (delta l_0)^2  
                        Kv=1.0_C_DOUBLE
                        denominator=denominator+visc_param*visc_param
                    case ( 33 ) !!Vatistas n=2 - visc_paramt<=>rc  
                        ! See Hoydonck 2012. and Script MainRingSensitivityN
                        ! Not matture enough
                        h = (norm_a+norm_b)/2.0_C_DOUBLE
                        visc_param2=visc_param**2
                        if(h<2*sqrt(norm2_r0)) then
                            ! use Vatistas n=2, normal one
                            h2 = norm2_crossprod/norm2_r0
                        end if
                        Kv = (h2)/sqrt(visc_param2**2+h2**2)
                    case DEFAULT
                        Kv=1.0_C_DOUBLE
                    end select 

                    Kv=Intensity*fourpi_inv*Kv*(norm_a+norm_b)/denominator
                    Uind(1:3) = Kv*crossprod(1:3)
                    !print "(A,3F21.16)","ui   :",Uind(1:3)

                    if (iComputeGrad==1) then 
                        Grad(1:9)=0.0_C_DOUBLE
                        !denom=(norm_a*norm_b + xa*xb+ya*yb+za*zb)+MINDENOM
                        !                     denom=denom+MINDENOM
                        !                     D=-(DP1/norm_a**3+DP2/norm_b**3)&
                        !                         -(1./norm_a+1./norm_b)*1/denom*(DP1/norm_a*norm_b+ DP2/norm_b*norm_a+ DP1+DP2)
                        ! 
                        !                     Grad(1) = crossprod(1)*D(1)                              ! 11
                        !                     Grad(2) = crossprod(1)*D(2) +(1./norm_a+1./norm_b)*(zb-za) ! 12
                        !                     Grad(3) = crossprod(1)*D(3) +(1./norm_a+1./norm_b)*(ya-yb) ! 13
                        ! 
                        !                     Grad(4) = crossprod(2)*D(1) +(1./norm_a+1./norm_b)*(za-zb) ! 21
                        !                     Grad(5) = crossprod(2)*D(2)                              ! 22
                        !                     Grad(6) = crossprod(2)*D(3) +(1./norm_a+1./norm_b)*(xb-xa) ! 23
                        ! 
                        !                     Grad(7) = crossprod(3)*D(1) +(1./norm_a+1./norm_b)*(yb-ya) ! 31
                        !                     Grad(8) = crossprod(3)*D(2) +(1./norm_a+1./norm_b)*(xa-xb) ! 32
                        !                     Grad(9) = crossprod(3)*D(3)                              ! 33
                        ! 
                        !                     Grad=Intensity*fourpi_inv*1./denom*Grad;
                        !                     if(visc_model>0) then
                        !                         ! h2 and visc_param2 are the same for everybody
                        !                         ! We use Lamb-Oseen model (it's the easiest)
                        !                         exp_value = -1.25643_C_DOUBLE*(h2)/(visc_param2)
                        !                         if(exp_value<MIN_EXP_VALUE) then
                        !                             ! Far away, we then use inviscid Gradient from above.
                        !                         else
                        !                             ! The Lamb-Oseen kernel (may be different than above Kv)
                        !                             Kv2 = 1.0_C_DOUBLE-exp(exp_value)
                        !                             ! Computing the Grad of Kv2
                        !                             D=0._C_DOUBLE
                        !                             ! Grad(raxrb) ra x rb 
                        !                             ! !!!!! Watcth out, need the segment extremity coordinates, but we don't have it
                        !                             ! It's the opposite formula hence, since xa-xb=x1-x2
                        !                             D(1)=(zb-za)*crossprod(2)+(ya-yb)*crossprod(3);
                        !                             D(2)=(za-zb)*crossprod(1)+(xb-xa)*crossprod(3);
                        !                             D(3)=(yb-ya)*crossprod(1)+(xa-xb)*crossprod(2);
                        !                             D=-2._C_DOUBLE*1.25643_C_DOUBLE/(norm2_r0 *visc_param2)*exp(exp_value)*D;
                        !                             !                         print*,'Kv',Kv
                        !                             !                         print*,'Grad',Grad
                        ! 
                        !                             ! Grad_viscous=Kv2*Grad + u_invisc * D
                        !                             Grad=Kv2*Grad ! Intensity already in Grad and Ui
                        ! Since Ui contains already Kv, we have to divide by it. Alternative: store Ui inviscid
                        !                            print*,'Kv',Kv,visc_param,exp_value,Kv2! this print is here to remind that c code has no grad
                        ! also a singulartiy can occur with the grid.
                        !                             Grad(1) = Grad(1) + Uind(1)*D(1)/Kv
                        !                             Grad(2) = Grad(2) + Uind(1)*D(2)/Kv
                        !                             Grad(3) = Grad(3) + Uind(1)*D(3)/Kv
                        ! 
                        !                             Grad(4) = Grad(4) + Uind(2)*D(1)/Kv
                        !                             !                         print*, 'Term3', Uind(2)*D(1)/Kv
                        !                             Grad(5) = Grad(5) + Uind(2)*D(2)/Kv
                        !                             Grad(6) = Grad(6) + Uind(2)*D(3)/Kv
                        ! 
                        !                             Grad(7) = Grad(7) + Uind(3)*D(1)/Kv
                        !                             Grad(8) = Grad(8) + Uind(3)*D(2)/Kv
                        !                             Grad(9) = Grad(9) + Uind(3)*D(3)/Kv
                        !                         endif
                        !                     end if  ! Visc model >0
                    else
                        Grad(1:9)=0.0_C_DOUBLE
                    end if  ! iComputeGrad
                end if ! denominator size or distances too small
            end if ! 
            ! --------------------------------------------------------------------------------
            ! --- END inlining of Segment 11 
            ! --------------------------------------------------------------------------------
            !icp-icp_beg+1
            Uind_out(1:3,ioff+icp-icp_beg) = Uind_out(1:3,ioff+icp-icp_beg)+Uind(1:3)
            if(iComputeGrad==1) then
                Grad_out(1:9,ioff+icp-icp_beg) = Grad_out(1:9,ioff+icp-icp_beg)+Grad(1:9)
            endif
        end do
    enddo
    !$OMP END DO 
    !$OMP END PARALLEL


end subroutine 



<|start_filename|>welib/airfoils/Makefile<|end_filename|>
all: test

test:
	python -m unittest -v


<|start_filename|>welib/vortilib/elements/fortran/MathConstants.f90<|end_filename|>
!> 
module MathConstants
    use PrecisionMod, only: MK, SP
    implicit none
    real(MK), parameter :: eps           = 2.220446049250313e-16_MK
    ! Pi related
    real(MK), parameter :: pi            = 3.141592653589793238462643383279502884197169399375105820974944_MK
    real(MK), parameter :: pi_half       = 1.570796326794896619231321691639751442098584699687552910487472_MK
    real(MK), parameter :: twopi         = 6.283185307179586476925286766559005768394338798750211641949889_MK
    real(MK), parameter :: fourpi        = 12.56637061435917295385057353311801153678867759750042328389977_MK
    real(MK), parameter :: twopi_inv     = 0.159154943091895335768883763372514362034459645740456448747667_MK
    real(MK), parameter :: fourpi_inv    = 0.079577471545947667884441881686257181017229822870228224373833_MK
    real(MK), parameter :: four_third_pi = 4.188790204786390984616857844372670512262892532500141094633259_MK
    real(MK), parameter :: rad2deg       = 57.29577951308232087679815481410517033240547246656432154916024_MK
    real(MK), parameter :: deg2rad       = 0.017453292519943295769236907684886127134428718885417254560971_MK
    ! NaN
    integer, parameter  :: NaN=-999999
    real(MK), parameter :: rNaN=-999999._MK
    real(SP), parameter :: NaN_SP=-999999._SP

    interface is_nan
        module procedure iis_nan, ris_nan
    end interface
contains
    logical function iis_nan(x)
        integer,intent(in) ::x
        iis_nan=(x==-NaN)
    end function
    
    logical function ris_nan(x)
        use PrecisionMod, only:PRECISION_EPS
        real(MK),intent(in) ::x
        ris_nan=abs(x-rNaN)<PRECISION_EPS
    end function


end module MathConstants



<|start_filename|>welib/vortilib/elements/fortran/Part.f90<|end_filename|>
!>  Takes nPart particle and nCPs control points
subroutine fui_part_c(CPs,ioff,iCP_beg, nCPs_loc, nCPs, iComputeGrad, Uind_out, Grad_out, &
        nPartTot, Part, iPartCurPos, PartSmoothModel, PartSmoothParam, PartIntensities) 
    use PrecisionMod, only: C_INT, C_DOUBLE
    implicit none
    ! Arguments declarations 
    integer(C_INT), intent(in)                        :: iComputeGrad    !< 
    integer(C_INT), intent(in)                        :: ioff
    integer(C_INT), intent(in)                        :: iCP_beg
    integer(C_INT), intent(in)                        :: nCPs_loc
    integer(C_INT), intent(in)                        :: nCPs
    real(C_DOUBLE), dimension(3,nCPs), intent(in)     :: CPs             !< 
    real(C_DOUBLE), dimension(3,nCPs), intent(inout)    :: Uind_out        !< Induced velocity vector - No Side effects!!!
    real(C_DOUBLE), dimension(9,nCPs), intent(inout)    :: Grad_out        !< Gradient vector - No Side effects
    integer(C_INT), intent(in)                        :: nPartTot
    real(C_DOUBLE), dimension(3,nPartTot), intent(in) :: Part
    integer(C_INT) ,intent(in)                        :: iPartCurPos
    integer(C_INT), intent(in)                        :: PartSmoothModel
    real(C_DOUBLE), dimension(nPartTot), intent(in)   :: PartSmoothParam
    real(C_DOUBLE), dimension(3,nPartTot), intent(in) :: PartIntensities
    integer :: icp,ip

    ! 
    real(C_DOUBLE),parameter :: MINNORM=1e-4_C_DOUBLE
    real(C_DOUBLE), parameter :: fourpi_inv    = 0.079577471545947667884441881686257181017229822870228224373833_C_DOUBLE
    ! --------------------------------------------------------------------------------
    ! ---  Variables for particle 11
    ! --------------------------------------------------------------------------------
    ! Kindof arguments
    real(C_DOUBLE), dimension(9) :: Grad       !< no side effects
    real(C_DOUBLE), dimension(3) :: Uind       !< no side effects
    real(C_DOUBLE), dimension(3) :: DeltaP     !< CP-PP "control point - particle point"
    real(C_DOUBLE), dimension(3) :: Omega      !< 
    real(C_DOUBLE)               :: Eps        !< 
    ! 
    real(C_DOUBLE),dimension(3)  :: C          !< Cross product of Omega and r
    real(C_DOUBLE)               :: E          !< Exponential poart for the mollifider
    real(C_DOUBLE)               :: Ebar       !< 1-E
    real(C_DOUBLE)               :: r3_inv     !< 
    real(C_DOUBLE)               :: n5_inv     !< 
    real(C_DOUBLE)               :: n2_inv     !< TODO, remove this one or r3..but watch out, look for the Grad trick of model 2
    real(C_DOUBLE)               :: rDeltaP    !< norm , distance between point and particle
    real(C_DOUBLE)               :: ScalarPart !< the part containing the inverse of the distance, but not 4pi, Mollifier

    
    !print*,'Part fortran',(iCP_beg+1),(iCP_beg+1)+nCPS_loc-1,nCPs
    !$OMP PARALLEL default(shared)
    !$OMP do private(&
    !$OMP& icp,ip,Uind,Grad,DeltaP,Omega,Eps,C,E,Ebar,r3_inv,n5_inv,n2_inv,rDeltaP,ScalarPart&
    !$OMP& ) schedule(runtime)
    ! loop on CPs 
    do icp=(iCP_beg+1),(iCP_beg+1)+nCPS_loc-1
        ! loop on particles
        do ip=1,iPartCurPos
            ! Setting up, kindof arguments for particle 11
            Uind(1:3) = 0.0_C_DOUBLE
            Grad(1:9) = 0.0_C_DOUBLE
            DeltaP = CPs(:,icp)-Part(1:3,ip)
            Eps    = PartSmoothParam(ip)
            Omega  = PartIntensities(1:3,ip)
            ! --------------------------------------------------------------------------------
            ! --- Inlining Particle 11 
            ! --------------------------------------------------------------------------------
            rDeltaP=sqrt(DeltaP(1)**2+ DeltaP(2)**2+ DeltaP(3)**2)! norm
            ! TODO, escape should depend on the 
            if (rDeltaP<MINNORM) then! to avoid singularity, save a bit of computation and use limit values
                !--------------------------------------------------------------------------------
                !--- Exactly on the Singularity 
                !--------------------------------------------------------------------------------
                !return

                ! NEW setting the gradient to 0, to avoid self influence (artifact of the low order Smoothing Kernel)
                ! => Commenting what is below
                !             if( SmoothModel==0) then 
                !                 !return
                !             end if
                !             if(iComputeGrad) then
                !                 if( SmoothModel==2) then 
                !                     r3_inv=1_C_DOUBLE/sqrt(Eps) ! Epsilon^3
                !                 else
                !                     r3_inv=1_C_DOUBLE/Eps ! Epsilon ^3
                !                 endif
                !                 ! For model 1 and 2 is the same limit
                !                 Grad(2) = fourpi_inv*( -Omega(3) *r3_inv) ! 12
                !                 Grad(3) = fourpi_inv*( +Omega(2) *r3_inv) ! 13
                !                 Grad(4) = fourpi_inv*( +Omega(3) *r3_inv) ! 21
                !                 Grad(6) = fourpi_inv*( -Omega(1) *r3_inv) ! 23
                !                 Grad(7) = fourpi_inv*( -Omega(2) *r3_inv) ! 31
                !                 Grad(8) = fourpi_inv*( +Omega(1) *r3_inv) ! 33
                !             endif
                ! TODO far distance
            else
                !--------------------------------------------------------------------------------
                !--- Normal Procedure 
                !--------------------------------------------------------------------------------


                C(1) = Omega(2) * DeltaP(3) - Omega(3) * DeltaP(2)
                C(2) = Omega(3) * DeltaP(1) - Omega(1) * DeltaP(3)
                C(3) = Omega(1) * DeltaP(2) - Omega(2) * DeltaP(1)
                ! influence of all control points 
                E    = 0.0_C_DOUBLE
                Ebar = 1.0_C_DOUBLE
                select case (PartSmoothModel) !
                case (0) ! No mollification, exept on zero
                    r3_inv=1._C_DOUBLE/(rDeltaP**3)
                    ScalarPart=r3_inv
                case (2) !
                    !Exponential mollifier
                    r3_inv=1._C_DOUBLE/(rDeltaP**3)
                    E=exp(-rDeltaP**3/Eps**3)
                    Ebar=1._C_DOUBLE-E
                    ScalarPart=Ebar*r3_inv
                case (3) !
                    !Exponential mollifier to the power 4
                    r3_inv=1._C_DOUBLE/(rDeltaP**3)
                    E=exp(-rDeltaP**4/Eps**4)
                    Ebar=1._C_DOUBLE-E
                    ScalarPart=Ebar*r3_inv
                case (22) !
                    ! Finite support
                    r3_inv= 1._C_DOUBLE/sqrt(Eps**6+rDeltaP**6)
                    ScalarPart= r3_inv
                case default 
                    print*,'SmoothModel',PartSmoothModel
                    print*,'Blob: wrong viscous model'
                    STOP
                end select 

                Uind(1:3)=C*ScalarPart*fourpi_inv

                !          write(*,*) 'DeltaP', DeltaP
                !          write(*,*) 'Omega', Omega
                !          write(*,*) 'Smoothparam', Eps
                !          write(*,*) 'SmoothModel', SmoothModel
                !          write(*,*) 'C', C
                !          write(*,*) 'r3_inv', r3_inv
                !          write(*,*) 'ScalarPart', ScalarPart
                !          write(*,*) 'UIout', UIout
                !          STOP
                if (iComputeGrad==1) then 
                    n5_inv=1._C_DOUBLE/rDeltaP**5
                    select case ( PartSmoothModel) !
                    ! Omega is my alpha
                    ! C = (alpha x r)  The Omega term comes from grad(alpha x r)
                    case (0) !
                        Grad(1) = fourpi_inv*(                   -3.*DeltaP(1) * C(1) *n5_inv) ! 11
                        Grad(2) = fourpi_inv*( -Omega(3) *r3_inv -3.*DeltaP(2) * C(1) *n5_inv) ! 12
                        Grad(3) = fourpi_inv*( +Omega(2) *r3_inv -3.*DeltaP(3) * C(1) *n5_inv) ! 13

                        Grad(4) = fourpi_inv*( +Omega(3) *r3_inv -3.*DeltaP(1) * C(2) *n5_inv) ! 21
                        Grad(5) = fourpi_inv*(                   -3.*DeltaP(2) * C(2) *n5_inv) ! 22
                        Grad(6) = fourpi_inv*( -Omega(1) *r3_inv -3.*DeltaP(3) * C(2) *n5_inv) ! 23

                        Grad(7) = fourpi_inv*( -Omega(2) *r3_inv -3.*DeltaP(1) * C(3) *n5_inv) ! 31
                        Grad(8) = fourpi_inv*( +Omega(1) *r3_inv -3.*DeltaP(2) * C(3) *n5_inv) ! 33
                        Grad(9) = fourpi_inv*(                   -3.*DeltaP(3) * C(3) *n5_inv) ! 32

                    case (2) !
                        ! Exponential mollifier 
                        n2_inv=1._C_DOUBLE/rDeltaP**2
                        ! Almost the same as above 
                        ! Remember I chose Eps=epsilon**3, so don't cube it twice!!
                        Grad(1) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 3./Eps**3*n2_inv*E)*DeltaP(1)*C(1))!11
                        Grad(2) = fourpi_inv*( -Omega(3)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./Eps**3*n2_inv*E)*DeltaP(2)*C(1))!12
                        Grad(3) = fourpi_inv*( +Omega(2)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./Eps**3*n2_inv*E)*DeltaP(3)*C(1))!13

                        Grad(4) = fourpi_inv*( +Omega(3)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./Eps**3*n2_inv*E)*DeltaP(1)*C(2))!21
                        Grad(5) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 3./Eps**3*n2_inv*E)*DeltaP(2)*C(2))!22
                        Grad(6) = fourpi_inv*( -Omega(1)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./Eps**3*n2_inv*E)*DeltaP(3)*C(2))!23

                        Grad(7) = fourpi_inv*( -Omega(2)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./Eps**3*n2_inv*E)*DeltaP(1)*C(3))!31
                        Grad(8) = fourpi_inv*( +Omega(1)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./Eps**3*n2_inv*E)*DeltaP(2)*C(3))!33
                        Grad(9) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 3./Eps**3*n2_inv*E)*DeltaP(3)*C(3))!32

                    case (3) !
                        ! Exponential mollifier 
                        n2_inv=1._C_DOUBLE/rDeltaP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!FAKE
                        ! Almost the same as above 
                        ! Remember I chose Eps=epsilon**3, so don't cube it twice!!
                        Grad(1) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 4./Eps**4*n2_inv*E)*DeltaP(1)*C(1))!11
                        Grad(2) = fourpi_inv*( -Omega(3)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./Eps**4*n2_inv*E)*DeltaP(2)*C(1))!12
                        Grad(3) = fourpi_inv*( +Omega(2)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./Eps**4*n2_inv*E)*DeltaP(3)*C(1))!13

                        Grad(4) = fourpi_inv*( +Omega(3)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./Eps**4*n2_inv*E)*DeltaP(1)*C(2))!21
                        Grad(5) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 4./Eps**4*n2_inv*E)*DeltaP(2)*C(2))!22
                        Grad(6) = fourpi_inv*( -Omega(1)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./Eps**4*n2_inv*E)*DeltaP(3)*C(2))!23

                        Grad(7) = fourpi_inv*( -Omega(2)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./Eps**4*n2_inv*E)*DeltaP(1)*C(3))!31
                        Grad(8) = fourpi_inv*( +Omega(1)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./Eps**4*n2_inv*E)*DeltaP(2)*C(3))!33
                        Grad(9) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 4./Eps**4*n2_inv*E)*DeltaP(3)*C(3))!32

                    case (22) !
                        ! Finite support with Epsilon 6 
                        n5_inv= r3_inv**3 * rDeltaP**4 ! Fake n5_inv
                        ! EXACTLY the same as case 0 due to fake variables n3 and n5 
                        Grad(1) = fourpi_inv*(                   -3.*DeltaP(1) * C(1) *n5_inv) ! 11
                        Grad(2) = fourpi_inv*( -Omega(3) *r3_inv -3.*DeltaP(2) * C(1) *n5_inv) ! 12
                        Grad(3) = fourpi_inv*( +Omega(2) *r3_inv -3.*DeltaP(3) * C(1) *n5_inv) ! 13

                        Grad(4) = fourpi_inv*( +Omega(3) *r3_inv -3.*DeltaP(1) * C(2) *n5_inv) ! 21
                        Grad(5) = fourpi_inv*(                   -3.*DeltaP(2) * C(2) *n5_inv) ! 22
                        Grad(6) = fourpi_inv*( -Omega(1) *r3_inv -3.*DeltaP(3) * C(2) *n5_inv) ! 23

                        Grad(7) = fourpi_inv*( -Omega(2) *r3_inv -3.*DeltaP(1) * C(3) *n5_inv) ! 31
                        Grad(8) = fourpi_inv*( +Omega(1) *r3_inv -3.*DeltaP(2) * C(3) *n5_inv) ! 33
                        Grad(9) = fourpi_inv*(                   -3.*DeltaP(3) * C(3) *n5_inv) ! 32
                    end select  ! Smooth Model
                else
                    Grad(1:9)=0.0_C_DOUBLE
                end if  ! iComputeGrad
                !print*,'Uind',Uind
            end if  ! Norm Big enough
            ! --------------------------------------------------------------------------------
            ! --- END of particle 11 inlining 
            ! --------------------------------------------------------------------------------
            !icp-icp_beg+1
            Uind_out(1:3,ioff+icp-icp_beg)=Uind_out(1:3,ioff+icp-icp_beg)+Uind(1:3)
            if(iComputeGrad==1) then
                Grad_out(1:9,ioff+icp-icp_beg)=Grad_out(1:9,ioff+icp-icp_beg)+Grad(1:9)
            endif
        enddo! loop particles
    enddo ! loop Cps
    !$OMP END DO 
    !$OMP END PARALLEL
end subroutine



<|start_filename|>welib/vortilib/elements/fortran/SurfaceMap.f90<|end_filename|>
!> 
module SurfaceMap
    use PrecisionMod, only: MK
    implicit none
            
contains
            
    ! function ue=fSurfaceIntegral(vX0,vY0,vZ0,Ux,Uy,Uz,n,XCP,YCP,ZCP)
    ! Compute the integral of the source and doublet from a given "Velocity on a Face" using point approximation
    subroutine fUi_PointSourceVortex_surface_mesh(P0,Sigmas,Gammas,Mesh,h)
        !use MathConstants, only: fourpi_inv
        use MeshTypes, only: T_Mesh
        use UIBlobs, only: fUi_Particle_11
        use UIPointSource, only: fUi_PointSource_11
        ! Arguments
        type(T_Mesh), intent(inout)         :: Mesh
        real(MK), dimension(:,:),intent(in) :: P0     !< Position of sources and doublet
        real(MK), dimension(:),intent(in)   :: Sigmas
        real(MK), dimension(:,:),intent(in) :: Gammas
        real(MK),intent(in)                 :: h
        ! Variables
        integer :: im1,im2,im3
        integer :: nm1,nm2,nm3
        integer :: p
        real(MK), dimension(3) :: CP
        real(MK), dimension(3) :: r
        real(MK), dimension(3) :: us
        real(MK), dimension(3) :: uv
        real(MK), dimension(3) :: Uind_cum
        real(MK), dimension(9) :: Grad_cum
        real(MK), dimension(3) :: Uind_tmp
        real(MK), dimension(9) :: Grad_tmp
        real(MK) :: r_norm
        real(MK) :: r3_inv
        real(MK) :: r2_inv
        real(MK) :: factor
        integer :: isd,nsd
        integer :: SmoothModel
        real(MK) :: SmoothParam

        ! Init
        nm1=Mesh%GP%n1
        nm2=Mesh%GP%n2
        nm3=Mesh%GP%n3
        nsd=size(P0,2)

        SmoothParam=1.5*h;

        !$OMP PARALLEL DEFAULT(SHARED)
        !$OMP DO PRIVATE(p,im1,im2,im3,CP,Uind_cum,Grad_cum,Uind_tmp,Grad_tmp,&
        !$OMP& isd,r,r_norm,r3_inv,r2_inv,factor,SmoothModel,us,uv&
        !$OMP& ) schedule(runtime)
        ! Loop on CPs is unrolled
        do p=1,Mesh%nCPs
            ! Mesh indexes
            im1=mod((p-1),nm1)+1
            im2=mod((p-1)/nm1,nm2)+1
            im3=(p-1)/(nm1*nm2)+1
            !
            CP(1:3)=(/Mesh%GP%v1(im1),Mesh%GP%v2(im2),Mesh%GP%v3(im3)/)
            Uind_cum(1:3)=0.0_MK
            Grad_cum(1:9)=0.0_MK

            ! The 'if' below was a safety to avoid the points next to the singularities
            !if(im1>1 .and. im2>1 .and. im3>1 .and. im1<nm1 .and. im2<nm2 .and. im3<nm3) then
                ! Loop on Sources, unrolled
                do isd=1,nsd
                    ! ue=-(un(i01,i02)*r+cross(squeeze(ut(i01,i02,1:3)),r))/(4*pi*r_norm^3)*dS;
                    r(1:3)      = CP(1:3)-P0(1:3,isd)
                    r_norm = sqrt(r(1)*r(1)+r(2)*r(2)+r(3)*r(3))

                    !
                    SmoothModel=1;
                    Uind_tmp(1:3)=0.0_MK
                    Grad_tmp(1:9)=0.0_MK
                    call fUi_PointSource_11(r,Sigmas(isd), SmoothModel, SmoothParam,.true.,Uind_tmp, Grad_tmp)
                    Uind_cum(1:3)= Uind_cum(1:3)+Uind_tmp(1:3)
                    Grad_cum(1:9)= Grad_cum(1:9)+Grad_tmp(1:9)
                    ! 
                    SmoothModel=2;
                    Uind_tmp(1:3)=0.0_MK
                    Grad_tmp(1:9)=0.0_MK
                    call fUi_Particle_11(r, Gammas(1:3,isd), SmoothModel, SmoothParam, .true., Uind_tmp, Grad_tmp)
                    Uind_cum(1:3)= Uind_cum(1:3)+Uind_tmp(1:3)
                    Grad_cum(1:9)= Grad_cum(1:9)+Grad_tmp(1:9)
!                     if (r_norm>h) then
!                         us(1) = Sigmas(isd)*r(1)  ! 1/(4*pi*r_norm^3)
!                         us(2) = Sigmas(isd)*r(2)  ! 1/(4*pi*r_norm^3)
!                         us(3) = Sigmas(isd)*r(3)  ! 1/(4*pi*r_norm^3)
!                         uv(1) = Gammas(2,isd) * r(3) - Gammas(3,isd) * r(2)
!                         uv(2) = Gammas(3,isd) * r(1) - Gammas(1,isd) * r(3)
!                         uv(3) = Gammas(1,isd) * r(2) - Gammas(2,isd) * r(1)
!                         r3_inv=1._MK/(r_norm**3);
!                         factor=fourpi_inv*r3_inv !1/(4 pi r^3)
!                         Uind_cum(1:3)= Uind_cum(1:3)+(us(1:3)+uv(1:3))*factor
! 
!                         r2_inv=1._MK/(r_norm**2);
!                         Grad_cum(1) = Grad_cum(1)+factor*Sigmas(isd)*(1 -3*r(1)*r(1)*r2_inv) ! 11
!                         Grad_cum(2) = Grad_cum(2)+factor*Sigmas(isd)*(  -3*r(2)*r(1)*r2_inv) ! 12
!                         Grad_cum(3) = Grad_cum(3)+factor*Sigmas(isd)*(  -3*r(3)*r(1)*r2_inv) ! 13
!                                                         
!                         Grad_cum(4) = Grad_cum(4)+factor*Sigmas(isd)*(  -3*r(1)*r(2)*r2_inv) ! 21
!                         Grad_cum(5) = Grad_cum(5)+factor*Sigmas(isd)*(1 -3*r(2)*r(2)*r2_inv) ! 22
!                         Grad_cum(6) = Grad_cum(6)+factor*Sigmas(isd)*(  -3*r(3)*r(2)*r2_inv) ! 23
!                                                         
!                         Grad_cum(7) = Grad_cum(7)+factor*Sigmas(isd)*(  -3*r(1)*r(3)*r2_inv) ! 31
!                         Grad_cum(8) = Grad_cum(8)+factor*Sigmas(isd)*(  -3*r(2)*r(3)*r2_inv) ! 32
!                         Grad_cum(9) = Grad_cum(9)+factor*Sigmas(isd)*(1 -3*r(3)*r(3)*r2_inv) ! 33
! 
!                     endif
                enddo
                Mesh%Values(1:3,im1,im2,im3)  = Mesh%Values(1:3,im1,im2,im3)+Uind_cum(1:3)
                Mesh%Values(4:12,im1,im2,im3) = Mesh%Values(4:12,im1,im2,im3)+Grad_cum(1:9)
            !endif
        enddo
        !$OMP END DO 
        !$OMP END PARALLEL
    end subroutine
end module SurfaceMap


<|start_filename|>welib/vortilib/elements/fortran/UIVortexRings.f90<|end_filename|>
!> 
module UIVortexRings
    use PrecisionMod, only: MK
    implicit none

    real(MK),parameter :: MINNORM=1e-4
contains

    !     subroutine fUi_VortexRings(CPs,PRingCenter,R,GammaOriented, SmoothModel,SmoothParam,bComputeGrad,Uind,Grad)
    !         use MathConstants, only: fourpi_inv
    !         ! Arguments declarations 
    !         integer, intent(in) :: nCPs
    !         integer, intent(in) :: nRings
    !         logical, intent(in) :: bComputeGrad !< 
    !         real(MK), dimension(3,nCPs), intent(in)    :: CPs         !< side effects
    !         real(MK), dimension(3,nCPs), intent(inout) :: Uind        !< side effects
    !         real(MK), dimension(9,nCPs), intent(inout) :: Grad        !< side effects
    !         real(MK), dimension(3,nRings), intent(in)  :: PRingCenter !< 
    !         real(MK), dimension(nSrc), intent(in)      :: Intensities !< 
    !         integer, intent(in)                        :: SmoothModel !< 
    !         real(MK), dimension(nSrc), intent(in)      :: SmoothParam !< 
    ! !         ! Variable declarations 
    !     end subroutine
    !
    subroutine fUi_VortexRing_11(CP,PRingCenter,R,GammaOriented, SmoothModel,SmoothParam,bComputeGrad,Uind,Grad)
        ! Arguments declarations 
        real(MK), dimension(3), intent(in)    :: CP           !< 
        real(MK), dimension(3), intent(in)    :: PRingCenter  !< PRingCenter : Position of ring center 
        real(MK), intent(in)                  :: R            !< Ring radius
        real(MK), dimension(3), intent(in)    :: GammaOriented!< Vector norm is ring intensity, direction is normal to the ring
        integer, intent(in)                   :: SmoothModel  !< 
        real(MK), intent(in)                  :: SmoothParam  !< 
        logical, intent(in)                   :: bComputeGrad !< 
        real(MK), dimension(3), intent(inout) :: Uind         !< side effects
        real(MK), dimension(9), intent(inout) :: Grad         !< side effects
        ! Variable declarations 
        real(MK), dimension(3) :: ex  !<  
        real(MK), dimension(3) :: ey  !< 
        real(MK), dimension(3) :: ez  !<  
        real(MK), dimension(3,3) :: T_r2g
        real(MK), dimension(3,3) :: T_g2r
        real(MK) :: Gamma  !<  
        integer :: imin
        real(MK), dimension(3) :: CP_loc  !< 
        real(MK) :: R_cp
        real(MK) :: Z_cp
        real(MK) :: phi
        real(MK), dimension(3) :: Uindtmp         !< 
        real(MK), dimension(9) :: Gradtmp         !<
        real(MK) :: urad
        real(MK) :: upsi
        !         ! Note: theories return only vr (index 1) and vz (index 3), there is some cartesian handling though, but does it work ? Is it consistent with the rest of the codes. Does it make use of the inputed vpsi to reduce the error in the angle estimate. 
        ! 
        !! Setting up an orthonornal cartesian system normal to the ring
        Gamma=norm2(GammaOriented)
        ez(1:3)=GammaOriented(1:3)/Gamma
        imin=minloc(abs(ez),1)
        ey=ez
        ey(imin)=0
        ey=ey([ 2,3,1 ] )
        ey(imin)=-ey(imin)
        ey=ey/norm2(ey)
        ex=cross(ey,ez)
        !! Now setting up rotation matrix 
        T_r2g(1:3,1)=ex(1:3)
        T_r2g(1:3,2)=ey(1:3)
        T_r2g(1:3,3)=ez(1:3)
        T_g2r=transpose(T_r2g)

        !! Coordinate transformation to a cartesian frame normal to the ring
        CP_loc(1)=T_g2r(1,1)*(CP(1)-PRingCenter(1)) + T_g2r(1,2)*(CP(2)-PRingCenter(2)) + T_g2r(1,3)*(CP(3)-PRingCenter(3))
        CP_loc(2)=T_g2r(2,1)*(CP(1)-PRingCenter(1)) + T_g2r(2,2)*(CP(2)-PRingCenter(2)) + T_g2r(2,3)*(CP(3)-PRingCenter(3))
        CP_loc(3)=T_g2r(3,1)*(CP(1)-PRingCenter(1)) + T_g2r(3,2)*(CP(2)-PRingCenter(2)) + T_g2r(3,3)*(CP(3)-PRingCenter(3))

        !! Coordinate transform to a canonical polar coordinates
        Z_cp=CP_loc(3) ! Since we are already in the ring coordinate system 
        R_cp=sqrt(CP_loc(1)**2+CP_loc(2)**2)
        !print*,'CP_loc',CP_loc

        !! Computing velocity in polar coordinates
        Uindtmp=0.0_MK
        Gradtmp=0.0_MK
        if(SmoothModel==0) then
            call fUi_Ring11_loc(R_cp, Z_cp, Gamma, R, Uindtmp)
            !print*,'Uind pol',Uindtmp
        else
            print*,'Smooth param',SmoothParam ! for unused var compiler warning
            write(*,*)'Smooth Ring not done'
            STOP -1
        endif

        !! Converting to cartesian 
        phi=atan2(CP_loc(2),CP_loc(1)) 
        !print*,'phi',phi*180/3.141593
        urad=Uindtmp(1)
        upsi=0.0_MK
        Uindtmp(1)=+urad*cos(phi)
        Uindtmp(2)=+urad*sin(phi)
        !print*,'Uind cart',Uindtmp

        !! Going back to global frame 
        Uind(1)=Uind(1) + T_r2g(1,1)*Uindtmp(1) + T_r2g(1,2)*Uindtmp(2) + T_r2g(1,3)*Uindtmp(3)
        Uind(2)=Uind(2) + T_r2g(2,1)*Uindtmp(1) + T_r2g(2,2)*Uindtmp(2) + T_r2g(2,3)*Uindtmp(3)
        Uind(3)=Uind(3) + T_r2g(3,1)*Uindtmp(1) + T_r2g(3,2)*Uindtmp(2) + T_r2g(3,3)*Uindtmp(3)

        !print*,'Uind glob',Uind
        !! TODO gradient 
        if (bComputeGrad) then 
            Grad = 0.0D0  
        else
            Grad = 0.0D0 
        end if 
    contains
        !> Cross product of two vectors ! See crossprod in UTILS...
        function cross(a, b)
            real(MK), dimension(3) :: cross
            real(MK), dimension(3), intent(in) :: a, b
            cross(1) = a(2) * b(3) - a(3) * b(2)
            cross(2) = a(3) * b(1) - a(1) * b(3)
            cross(3) = a(1) * b(2) - a(2) * b(1)
        end function cross
    end subroutine

    subroutine test_ring()
        ! Arguments declarations 
        real(MK):: Gamma !< Ring intensity
        real(MK):: R     !< Ring radius
        real(MK), dimension(3) :: Uind  !< side effects, in coord (r, theta, z)
        real(MK), dimension(3) :: CP           !< 
        real(MK), dimension(3):: PRingCenter  !< PRingCenter : Position of ring center 
        real(MK), dimension(3):: GammaOriented!< Vector norm is ring intensity, direction is normal to the ring
        integer:: SmoothModel  !< 
        real(MK):: SmoothParam  !< 
        logical:: bComputeGrad !< 
        real(MK), dimension(9):: Grad         !< side effects
       
        R=3._MK
        Gamma=2._MK
        Uind=0_MK
        call fUi_Ring11_loc(1._MK, 2._MK, Gamma, R, Uind)
        GammaOriented=[2,0,0];
        PRingCenter=[3,2,1];
        CP=[5,2,1]
        CP=CP+[0.,cos(0.1),sin(0.1)];

        SmoothModel=0
        SmoothParam=0.0
        bComputeGrad=.false.
        Uind=[0,0,0]
        Grad(1:9)=0
        call fUi_VortexRing_11(CP,PRingCenter,R,GammaOriented, SmoothModel,SmoothParam,bComputeGrad,Uind,Grad)
    end subroutine


    !> Induction from a ring at a control point, within the referential of a ring: coordinates (r,theta, z), the ring is at z=0
    subroutine fUi_Ring11_loc(Rcp, Zcp, Gamma, R, Uind)
        use EllipticIntegrals, only: elliptic_ke
        use MathConstants, only: fourpi_inv
        ! Arguments declarations 
        real(MK), intent(in) :: Rcp   !< Control point radius
        real(MK), intent(in) :: Zcp   !< Control point z location
        real(MK), intent(in) :: Gamma !< Ring intensity
        real(MK), intent(in) :: R     !< Ring radius
        real(MK), dimension(3), intent(inout) :: Uind  !< side effects, in coord (r, theta, z)
        ! Variable declarations 
        real(MK), parameter :: MIN_DIST=1e-7 ! threshold for using axis formula
        real(MK) :: a  !<  
        real(MK) :: AA !<  
        real(MK) :: B  !<  
        real(MK) :: E  !<  
        real(MK) :: I1  !< 
        real(MK) :: I2  !<
        real(MK) :: K  !<  
        real(MK) :: m  !<  
        !real(MK) :: aL,x,rr,cons,phi,a2  !<  Lewis notations

        if (Rcp<MIN_DIST*R) then
            !! Using Axis formula : v_z=-Gamma/(2R) *1 / (1+(z/R)^2)^(3/2) See Naca 1184 - Watsle, or Glauert Contraction of slipstream of an airscrew 1926 
            Uind(1)=0
            Uind(3)=Gamma/(2_MK*R)*(1_MK/((1_MK+(Zcp/R)**2)**(1.5_MK))) ! no minus sign... it's a matter of convention
        else
            ! Formulation uses Formula from Yoon 2004.. it's not my favorite though.. 
            !Rcp=sqrt(Xcp.**2+Ycp.**2)
            !Zcp=Zcp-z ! !!!! In this case, rememeber to remove z to be in Ring coordinate system
            a=sqrt((Rcp+R)**2+(Zcp)**2)
            m=4_MK*Rcp*R/(a**2)
            AA=(Zcp)**2+Rcp**2+R**2
            B=-2_MK*Rcp*R
            call elliptic_ke(m, K, E) 
            I1=4_MK/a*K
            I2=4_MK/a**3*E/(1_MK-m)
            Uind(1)=( (Zcp)/B )*(I1-AA*I2)      ! vr  !! Watchout GammaR/4pi removed , norm2(r_i) for yoon is R
            Uind(3)=((R+Rcp*AA/B)*I2-Rcp/B*I1 ) ! vz 
            ! General scaling 
            Uind(1:3)=Uind(1:3)*Gamma*fourpi_inv*R

            ! Lewis
!             x =Zcp/R;
!             rr=Rcp/R;
!             aL    = x**2+(rr-1._MK)**2
!             a2    = x**2+(rr+1._MK)**2
!             cons = 0.5_MK/(pi*R*sqrt(x**2+(rr+1._MK)**2))*Gamma
!             m    = 4._MK*rr/a2
!             phi=atan(sqrt(4._MK*rr/aL))
!             if(m>0.999) then
!                 K=log(4./cos(phi))
!                 E=1.+0.5*(K-1.0/1.2) * cos(phi)**2
!             else
!                 call elliptic_ke(m, K, E) 
!             endif
!             Uind(3)= -cons*     (K-(1.+2.*(rr-1.)/aL)*E);
!             Uind(1)=  cons*x/rr*(K-(1.+2.*     rr/aL)*E);
            ! Oye with Grad
!                 zc=Zcp_prime(i)
!                 rc=Rcp_prime(i)
!                 a=sqrt((rc+R)**2+(zc)**2)
!                 m=4*rc*R/a**2
!                 A=(zc)**2+rc**2+R**2
!                 B=-2*rc*R
!                 call ellipke( m , K,E) 
!                 I1=4/a*K
!                 I2=4/a**3*E/(1-m)
!                 uip(i,2)=0
!                 uip(i,1)=Gamma/(4*pi)*R*( (zc)/B )*(I1-A*I2) ! vr hopefully
!                 uip(i,3)=Gamma/(4*pi)*R*((R+rc*A/B)*I2-rc/B*I1 ) ! vz hopefully
!                 if (bComputeGrad) then 
!                     a=(R + rc)**2 + zc**2
!                     b=(R + rc)**2 + zc**2
!                     ! X and Z are the same 
!                     ! Y and R are the same 
!                     Uzz=(1/2).*a.**(-3/2).*b.**(-2).*pi.**(-1).*zc.*(((-7).*R.**4+6.*R.**2.*rc.**2+rc.**4+2.*((-3).*R.**2+rc.**2).*zc.**2+zc.**4).*E+(-1).*b.*((-1).*R.**2+rc.**2+zc.**2).*K)
! 
!                     Uzr=(1/2).*a.**(-3/2).*b.**(-2).*pi.**(-1).*rc.**(-1).*(((R.**2+(-1).*rc.**2).**2.*(R.**2+rc.**2)+2.*(R.**4+(-6).*R.**2.*rc.**2+rc.**4).*zc.**2+(R.**2+rc.**2).*zc.**4).*E+(-1).*b.*((R.**2+(-1).*rc.**2).**2+(R.**2+rc.**2).*zc.**2).*K)
! 
!                     Urr= (-1/2).*a.**(-3/2).*b.**(-2).*pi.**(-1).*rc.**(-2).*zc.*((2.*rc.**6+5.* rc.**4.*(R.**2+zc.**2)+4.*rc.**2.*((-2).*R.**2+zc.**2).*(R.**2+zc.**2)+( R.**2+zc.**2).**3).*E+(-1).*b.*(2.*rc.**4+3.*rc.**2.*((-1).* R+zc).*(R+zc)+(R.**2+zc.**2).**2).*K)
! 
!                     Urz=(1/2).*a.**(-3/2).*b.**(-2).*pi.**(-1).*rc.**(-1).*(((R.**2+(-1).* rc.**2).**2.*(R.**2+rc.**2)+2.*(R.**4+(-6).*R.**2.*rc.**2+rc.**4).*zc.**2+( R.**2+rc.**2).*zc.**4).*E+(-1).*b.*((R.**2+(-1).*rc.**2).**2+ (R.**2+rc.**2).*zc.**2).*K)
! 
!                     grad(i,1)=Urr*Gamma
!                     grad(i,3)=Urz*Gamma
!                     grad(i,7)=Uzr*Gamma
!                     grad(i,9)=Uzz*Gamma
!                 end if 




        endif
    end subroutine fUi_Ring11_loc

    subroutine fui_ring11_loc_c(Rcp, Zcp, Gamma, R, Uind) bind(C,name='fui_ring11_loc_c')
        use PrecisionMod, only: C_DOUBLE
        !DEC$ ATTRIBUTES DLLEXPORT :: fui_ring11_loc_c
        !GCC$ ATTRIBUTES DLLEXPORT :: fui_ring11_loc_c
        ! Arguments declarations 
        real(C_DOUBLE), intent(in) :: Rcp   !< Control point radius
        real(C_DOUBLE), intent(in) :: Zcp   !< Control point z location
        real(C_DOUBLE), intent(in) :: Gamma !< Ring intensity
        real(C_DOUBLE), intent(in) :: R     !< Ring radius
        real(C_DOUBLE), dimension(3), intent(inout) :: Uind  !< side effects, in coord (r, theta, z)
        real(MK), dimension(3) :: Uind_MK 

        Uind_MK=0.0_MK
        call fUi_Ring11_loc(Rcp, Zcp, Gamma, R, Uind_MK)
        Uind=real(Uind_MK,C_DOUBLE)

    end subroutine 

end module


<|start_filename|>welib/vortilib/elements/fortran/UIHillsVortex.f90<|end_filename|>
!> 
module HillsVortex
    use  PrecisionMod, only: MK
    implicit none
            
contains
            
    !> Velocity induced by Hills vortex (hence in the frame with no free stream at infinity). 
    ! Induced velocity returned in Cartesian coordinates
    !
    ! Cylindrical coordinate system (rho,phi,z)
    subroutine fUi_HillsVortex_1(CP,a,us,z0,Uind,Grad,Defm,Vort,Strm)
        ! Arguments
        real(MK),dimension(3), intent(in) :: CP !< Control Point
        real(MK), intent(in) :: a  !< hill vortex radius
        real(MK), intent(in) :: us !< hills self induced velocity us=-5/2 U0  (ie convects along -ez)
        real(MK), intent(in) :: z0 !< Position on the z axis
        real(MK),dimension(3),intent(out) :: Uind !< Induced velocity
        real(MK),dimension(9),intent(out) :: Grad !< Gradient
        real(MK),dimension(3),intent(out) :: Defm !< Deformation
        real(MK),dimension(3),intent(out) :: Vort !< Vorticity
        real(MK),dimension(3),intent(out),optional :: Strm !< Stremfunction (optional)
        ! Variables
        real(MK)              :: r3d,rho
        real(MK)              :: x,y,z
        real(MK)              :: urho
        real(MK)              :: uz
        real(MK),dimension(3) :: e_phi  ! Tangential coordinate
        real(MK),dimension(3) :: e_rho  ! Cylindrical radial
        real(MK)              :: om_phi
        real(MK)              :: defm_phi
        real(MK)              :: U0 !< I don't want to think : U0 =-2/5 us
        real(MK)              :: Stokes 
        real(MK)              :: Stokes0
        !
        U0=-2._MK/5._MK*us
        !
        !
        x=CP(1)
        y=CP(2)
        z=CP(3)-z0  ! Note the allowed offset

        r3d=sqrt(x**2 + y**2 + z**2)
        rho=sqrt(x**2 + y**2)
        ! Tangential coordinate (tangent to (y,z))
        if (rho/a<1.e-12) then
            ! This should save us from singularities
            e_phi(1:3)=0.0_MK
            e_rho(1:3)=0.0_MK
        else
            e_phi(1)=-y/rho
            e_phi(2)= x/rho
            e_phi(3)=0
            e_rho(1)= x/rho
            e_rho(2)= y/rho
            e_rho(3)=0
        endif
        if(r3d<a) then 
            ! --- Inside the sphere 
            ! Velocity
            uz     =3._MK/5._MK*us*(1._MK-(2._MK*rho**2+z**2)/(a**2))+us*2._MK/5._MK ;
            urho   =3._MK/5._MK*us*(rho*z)/(a**2);
            ! Vorticity (along e_phi)
            om_phi= 3._MK*us*rho/a**2  != -15/2  * (u0 rho)/a^2 = 3 us rho/a**2
            ! Deformation
            defm_phi= 9._MK/5._MK*us**2/a**2*rho*z ! =45/4 * uo^2/a^2 * rho z
            ! Gradient
            Grad(1:9)=0.0_MK !TODO
            ! Stokes Stream function
            Stokes=-3.0_MK/4.0_MK*U0*rho**2*(1-(rho**2+z**2)/(a**2))
            Stokes0=-rho**2*U0/2.0_MK !< Due to the addition of a free stream term in uz_in
        else 
            ! --- Outside of the sphere 
            ! Velocity
            uz    =2._MK/5._MK*us* (((a**2)/(z**2+rho**2))**(5._MK/2._MK))*(2._MK*z**2-rho**2)/(2._MK*a**2);
            urho  =3._MK/5._MK*us*rho*z/(a**2)*(((a**2)/(z**2+rho**2))**(5._MK/2._MK));
            ! Vorticity
            om_phi   = 0.0_MK
            ! Deformation
            defm_phi = 0.0_MK
            ! Gradient 
            Grad(1:9)=0.0_MK !TODO
            ! Stokes Streamfunction
            Stokes=1.0_MK/2.0_MK*U0*rho**2*(1-a**3/((rho**2+z**2)**1.5_MK))
            Stokes0=-rho**2*U0/2.0_MK !< Due to the addition of a free stream term in uz_in
        endif
        !
        Uind(1) = urho * e_rho(1)
        Uind(2) = urho * e_rho(2)
        Uind(3) = uz
        !
        Defm(1) =  defm_phi * e_phi(1)
        Defm(2) =  defm_phi * e_phi(2)
        Defm(3) = 0.0_MK
        !
        Vort(1) = om_phi * e_phi(1)
        Vort(2) = om_phi * e_phi(2)
        Vort(3) = 0.0_MK
        ! Streamfunction
        if(present(Strm)) then
            if (rho/a>1.e-12) then
                Strm(1) = (Stokes0+Stokes)/rho * e_phi(1)
                Strm(2) = (Stokes0+Stokes)/rho * e_phi(2)
                Strm(3) = 0.0_MK
            else
                Strm=0.0_MK
            endif
        endif

    end subroutine 

end module HillsVortex


<|start_filename|>welib/airfoils/matlab/Makefile<|end_filename|>
all: prog test

prog:
	gfortran KarmanTrefftzGeometry.f90 -o KarmanTrefftzGeometry

test:
	# xc, yc, tau_deg
	echo "-0.1 0. 5" | ./KarmanTrefftzGeometry
	mv geom*.dat data/


<|start_filename|>welib/vortilib/elements/fortran/UIBlobs.f90<|end_filename|>
!> Module for blobs induced velocities
module UIBlobs
    use PrecisionMod, only: MK
    implicit none

    real(MK),parameter :: MINNORM=1e-4

contains

    !>  Takes nPart particle and nCPs control points
    subroutine fUi_Particle(CPs,Part, Omega, SmoothModel, SmoothParam, bComputeGrad, Uiout, Gradout,nCPs,nPart)
        ! Arguments declarations 
        integer, intent(in) :: nCPs
        integer, intent(in) :: nPart
        logical, intent(in) :: bComputeGrad !< 
        real(MK), dimension(3,nCPs), intent(in) :: CPs   !< side effects
        real(MK), dimension(3,nCPs), intent(inout) :: UiOut   !< side effects
        real(MK), dimension(9,nCPs), intent(inout)  :: Gradout !< side effects
        real(MK), dimension(3,nPart), intent(in)   :: Part      !< 
        real(MK), dimension(3,nPart), intent(in)   :: Omega      !< 
        integer, intent(in)    :: SmoothModel !< 
        real(MK), dimension(nPart), intent(in)   :: SmoothParam !< 
        ! Variable declarations 
        real(MK), dimension(3) :: UItmp   !< 
        real(MK), dimension(9) :: Gradtmp !< 
        real(MK), dimension(3) :: DP      !< 
        integer :: icp,ip

        ! loop on CPs 
        do icp=1,nCPs
            ! loop on particles
            do ip=1,nPart
                UItmp=0.0_MK
                Gradtmp=0.0_MK 
                DP=CPs(:,icp)-Part(1:3,ip)
                call  fUi_Particle_11(DP, Omega(1:3,ip), SmoothModel ,&
                    SmoothParam(ip), bComputeGrad, UItmp,Gradtmp)
                UIout(1:3,icp)=UIout(1:3,icp)+UItmp(1:3)
                if(bComputeGrad) then
                    Gradout(1:9,icp)=Gradout(1:9,icp)+Gradtmp(1:9)
                endif
            enddo! loop particles
        enddo ! loop Cps
    end subroutine

    !>  Takes one particle and one point
    subroutine fUi_Particle_11(DeltaP, Omega, SmoothModel, SmoothParam, bComputeGrad, Ui, Grad)
        use MathConstants, only: fourpi_inv
        ! Arguments declarations 
        logical, intent(in) :: bComputeGrad  !<  
        real(MK), dimension(9), intent(out) :: Grad  !< no side effects
        real(MK), dimension(3), intent(out) :: Ui  !<  no side effects
        real(MK), dimension(3), intent(in) :: DeltaP  !<  CP-PP "control point - particle point"
        real(MK), dimension(3), intent(in) :: Omega  !<  
        integer, intent(in) :: SmoothModel  !<  
        real(MK), intent(in) :: SmoothParam  !<  
        ! Variable declarations 
        real(MK),dimension(3) :: C  !<  Cross product of Omega and r
        real(MK) :: E  !<   Exponential poart for the mollifider
        real(MK) :: Ebar  !<  1-E
        real(MK) :: r3_inv  !<  
        real(MK) :: n5_inv  !<  
        real(MK) :: n2_inv  !<   TODO, remove this one or r3..but watch out, look for the Grad trick of model 2
        real(MK) :: rDeltaP  !< norm , distance between point and particle 
        real(MK) :: ScalarPart  !< the part containing the inverse of the distance, but not 4pi, Mollifier

        Ui(1:3)   = 0.0_MK
        Grad(1:9) = 0.0_MK
        rDeltaP=sqrt(DeltaP(1)**2+ DeltaP(2)**2+ DeltaP(3)**2)! norm
        ! TODO, escape should depend on the 
        if (rDeltaP<MINNORM) then! to avoid singularity, save a bit of computation and use limit values
            !--------------------------------------------------------------------------------
            !--- Exactly on the Singularity 
            !--------------------------------------------------------------------------------
            return
        
            ! NEW setting the gradient to 0, to avoid self influence (artifact of the low order Smoothing Kernel)
            ! => Commenting what is below
            !             if( SmoothModel==0) then 
            !                 return
            !             end if
            !             if(bComputeGrad) then
            !                 if( SmoothModel==2) then 
            !                     r3_inv=1_MK/sqrt(SmoothParam) ! Epsilon^3
            !                 else
            !                     r3_inv=1_MK/SmoothParam ! Epsilon ^3
            !                 endif
            !                 ! For model 1 and 2 is the same limit
            !                 Grad(2) = fourpi_inv*( -Omega(3) *r3_inv) ! 12
            !                 Grad(3) = fourpi_inv*( +Omega(2) *r3_inv) ! 13
            !                 Grad(4) = fourpi_inv*( +Omega(3) *r3_inv) ! 21
            !                 Grad(6) = fourpi_inv*( -Omega(1) *r3_inv) ! 23
            !                 Grad(7) = fourpi_inv*( -Omega(2) *r3_inv) ! 31
            !                 Grad(8) = fourpi_inv*( +Omega(1) *r3_inv) ! 33
            !             endif
            ! TODO far distance
        else
            !--------------------------------------------------------------------------------
            !--- Normal Procedure 
            !--------------------------------------------------------------------------------


            C(1) = Omega(2) * DeltaP(3) - Omega(3) * DeltaP(2)
            C(2) = Omega(3) * DeltaP(1) - Omega(1) * DeltaP(3)
            C(3) = Omega(1) * DeltaP(2) - Omega(2) * DeltaP(1)
            ! influence of all control points 
            select case (SmoothModel) !
            case (0) ! No mollification, exept on zero
                r3_inv=1._MK/(rDeltaP**3)
                ScalarPart=r3_inv
            case (2) !
                !Exponential mollifier
                r3_inv=1._MK/(rDeltaP**3)
                E=exp(-rDeltaP**3/SmoothParam**3)
                Ebar=1._MK-E
                ScalarPart=Ebar*r3_inv
            case (3) !
                !Exponential mollifier to the power 4
                r3_inv=1._MK/(rDeltaP**3)
                E=exp(-rDeltaP**4/SmoothParam**4)
                Ebar=1._MK-E
                ScalarPart=Ebar*r3_inv
            case (22) !
                ! Finite support
                r3_inv= 1._MK/sqrt(SmoothParam**6+rDeltaP**6)
                ScalarPart= r3_inv
            case default 
                print*,'SmoothModel',SmoothModel
                print*,'Blob: wrong viscous model'
            end select 

            Ui(1:3)=C*ScalarPart*fourpi_inv

        !          write(*,*) 'DeltaP', DeltaP
        !          write(*,*) 'Omega', Omega
        !          write(*,*) 'Smoothparam', SmoothParam
        !          write(*,*) 'SmoothModel', SmoothModel
        !          write(*,*) 'C', C
        !          write(*,*) 'r3_inv', r3_inv
        !          write(*,*) 'ScalarPart', ScalarPart
        !          write(*,*) 'UIout', UIout
        !          STOP
            if (bComputeGrad) then 
                n5_inv=1._MK/rDeltaP**5
                select case ( SmoothModel) !
                ! Omega is my alpha
                ! C = (alpha x r)  The Omega term comes from grad(alpha x r)
                case (0) !
                    Grad(1) = fourpi_inv*(                   -3.*DeltaP(1) * C(1) *n5_inv) ! 11
                    Grad(2) = fourpi_inv*( -Omega(3) *r3_inv -3.*DeltaP(2) * C(1) *n5_inv) ! 12
                    Grad(3) = fourpi_inv*( +Omega(2) *r3_inv -3.*DeltaP(3) * C(1) *n5_inv) ! 13

                    Grad(4) = fourpi_inv*( +Omega(3) *r3_inv -3.*DeltaP(1) * C(2) *n5_inv) ! 21
                    Grad(5) = fourpi_inv*(                   -3.*DeltaP(2) * C(2) *n5_inv) ! 22
                    Grad(6) = fourpi_inv*( -Omega(1) *r3_inv -3.*DeltaP(3) * C(2) *n5_inv) ! 23

                    Grad(7) = fourpi_inv*( -Omega(2) *r3_inv -3.*DeltaP(1) * C(3) *n5_inv) ! 31
                    Grad(8) = fourpi_inv*( +Omega(1) *r3_inv -3.*DeltaP(2) * C(3) *n5_inv) ! 33
                    Grad(9) = fourpi_inv*(                   -3.*DeltaP(3) * C(3) *n5_inv) ! 32

                case (2) !
                    ! Exponential mollifier 
                    n2_inv=1._MK/rDeltaP**2
                    ! Almost the same as above 
                    ! Remember I chose SmoothParam=epsilon**3, so don't cube it twice!!
                    Grad(1) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 3./SmoothParam**3*n2_inv*E)*DeltaP(1)*C(1))!11
                    Grad(2) = fourpi_inv*( -Omega(3)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./SmoothParam**3*n2_inv*E)*DeltaP(2)*C(1))!12
                    Grad(3) = fourpi_inv*( +Omega(2)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./SmoothParam**3*n2_inv*E)*DeltaP(3)*C(1))!13

                    Grad(4) = fourpi_inv*( +Omega(3)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./SmoothParam**3*n2_inv*E)*DeltaP(1)*C(2))!21
                    Grad(5) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 3./SmoothParam**3*n2_inv*E)*DeltaP(2)*C(2))!22
                    Grad(6) = fourpi_inv*( -Omega(1)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./SmoothParam**3*n2_inv*E)*DeltaP(3)*C(2))!23

                    Grad(7) = fourpi_inv*( -Omega(2)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./SmoothParam**3*n2_inv*E)*DeltaP(1)*C(3))!31
                    Grad(8) = fourpi_inv*( +Omega(1)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 3./SmoothParam**3*n2_inv*E)*DeltaP(2)*C(3))!33
                    Grad(9) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 3./SmoothParam**3*n2_inv*E)*DeltaP(3)*C(3))!32

                case (3) !
                    ! Exponential mollifier 
                    n2_inv=1._MK/rDeltaP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!FAKE
                    ! Almost the same as above 
                    ! Remember I chose SmoothParam=epsilon**3, so don't cube it twice!!
                    Grad(1) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 4./SmoothParam**4*n2_inv*E)*DeltaP(1)*C(1))!11
                    Grad(2) = fourpi_inv*( -Omega(3)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./SmoothParam**4*n2_inv*E)*DeltaP(2)*C(1))!12
                    Grad(3) = fourpi_inv*( +Omega(2)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./SmoothParam**4*n2_inv*E)*DeltaP(3)*C(1))!13

                    Grad(4) = fourpi_inv*( +Omega(3)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./SmoothParam**4*n2_inv*E)*DeltaP(1)*C(2))!21
                    Grad(5) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 4./SmoothParam**4*n2_inv*E)*DeltaP(2)*C(2))!22
                    Grad(6) = fourpi_inv*( -Omega(1)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./SmoothParam**4*n2_inv*E)*DeltaP(3)*C(2))!23

                    Grad(7) = fourpi_inv*( -Omega(2)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./SmoothParam**4*n2_inv*E)*DeltaP(1)*C(3))!31
                    Grad(8) = fourpi_inv*( +Omega(1)*r3_inv*Ebar + (-3.*n5_inv*Ebar + 4./SmoothParam**4*n2_inv*E)*DeltaP(2)*C(3))!33
                    Grad(9) = fourpi_inv*(                       + (-3.*n5_inv*Ebar + 4./SmoothParam**4*n2_inv*E)*DeltaP(3)*C(3))!32

                case (22) !
                    ! Finite support with Epsilon 6 
                    n5_inv= r3_inv**3 * rDeltaP**4 ! Fake n5_inv
                    ! EXACTLY the same as case 0 due to fake variables n3 and n5 
                    Grad(1) = fourpi_inv*(                   -3.*DeltaP(1) * C(1) *n5_inv) ! 11
                    Grad(2) = fourpi_inv*( -Omega(3) *r3_inv -3.*DeltaP(2) * C(1) *n5_inv) ! 12
                    Grad(3) = fourpi_inv*( +Omega(2) *r3_inv -3.*DeltaP(3) * C(1) *n5_inv) ! 13

                    Grad(4) = fourpi_inv*( +Omega(3) *r3_inv -3.*DeltaP(1) * C(2) *n5_inv) ! 21
                    Grad(5) = fourpi_inv*(                   -3.*DeltaP(2) * C(2) *n5_inv) ! 22
                    Grad(6) = fourpi_inv*( -Omega(1) *r3_inv -3.*DeltaP(3) * C(2) *n5_inv) ! 23

                    Grad(7) = fourpi_inv*( -Omega(2) *r3_inv -3.*DeltaP(1) * C(3) *n5_inv) ! 31
                    Grad(8) = fourpi_inv*( +Omega(1) *r3_inv -3.*DeltaP(2) * C(3) *n5_inv) ! 33
                    Grad(9) = fourpi_inv*(                   -3.*DeltaP(3) * C(3) *n5_inv) ! 32
                end select  ! Smooth Model
            else
                Grad(1:9)=0.0_MK
            end if  ! bComputeGrad
        end if  ! Norm Big enough
    end subroutine fUi_Particle_11
end module UIBlobs



<|start_filename|>welib/vortilib/elements/fortran/UIVortexPoint2D.f90<|end_filename|>
!> 
module VortexPoint2D
    use PrecisionMod,   only: MK
    implicit none
    real(MK),parameter :: MINNORM2=1e-8
contains
    subroutine fUi_VortexPoint2D_11(DeltaP,Intensity,SmoothModel,SmoothParam,bComputeGrad,Uind, Grad)
        use MathConstants , only: twopi_inv
        implicit none 
        ! Input/output arguments 
        real(MK), dimension(3), intent(in) :: DeltaP       !< 3 x 1   Pcp-Pv
        real(MK), intent(in)               :: Intensity    !<  Circulation Gamma
        integer, intent(in) ,optional      :: SmoothModel  !< 
        real(MK), intent(in),optional      :: SmoothParam  !< 
        logical, intent(in)                :: bComputeGrad !< 
        real(MK), dimension(3),intent(out) :: Uind         !< No Side effects. Note that
        real(MK), dimension(9),intent(out) :: Grad         !< No Side effects
        !         ! Variables declaration 
        real(MK) :: r2     !< 
        real(MK) :: factor !< 
        real(MK) :: E      !< 
        real(MK) :: rho2   !< 
        real(MK) :: Qm     !< 
        ! 
        ! No need of gradient in 2D
        if(bComputeGrad) then
            Grad(1:9)=0.0_MK !
        else
            Grad(1:9)=0.0_MK
        endif
        Uind=0.0_MK

        r2  = DeltaP(1)**2+ DeltaP(2)**2
        if (r2<MINNORM2) then! to avoid strong singularity
            Uind(3)=0.0_MK
            Grad(1:9)=0.0_MK 
!         elseif (r2>100) then 
!             Uind(3)=0.0_MK
!             Grad(1:9)=0.0_MK 
        else
            select case (SmoothModel) !
            case (0) ! No mollification, exept on zero
                factor=1._MK/r2
            case (22) !
                ! Finite support 
                factor=1._MK/sqrt(SmoothParam**4+r2**2)
                ! --------------------------------------------------------------------------------
                ! --- Majda Kernels
                ! --------------------------------------------------------------------------------
            case (2) ! Second order - Gaussian - Majada
                rho2   = r2/SmoothParam**2;
                E      = exp(-rho2);
                Qm     = 1._MK;
                factor = (1._MK-Qm*E)/r2;
            case(4)  ! Fourth Order - Gaussian - Majda
                rho2   = r2/SmoothParam**2;
                E      = exp(-rho2);
                Qm     = 1._MK-rho2;
                factor = (1._MK-Qm*E)/r2;
            case(6)  ! Sixth Order - Gaussian - Majda
                rho2   = r2/SmoothParam**2;
                E      = exp(-rho2);
                Qm     = (1._MK-2._MK*rho2+rho2**2/2._MK);
                factor = (1._MK-Qm*E)/r2;
            case(8)  ! Eigth Order - Gaussian - Majda
                rho2   = r2/SmoothParam**2;
                E      = exp(-rho2);
                Qm     = (1._MK-3_MK*rho2+3._MK*rho2**2/2._MK-rho2**3/6._MK);
                factor = (1._MK-Qm*E)/r2;
            case(10)  ! Tenth Order - Gaussian - Mads
                rho2   = r2/SmoothParam**2;
                E      = exp(-rho2/2._MK);  !!!!!!!! Note the 2 in the Gaussian
                Qm     = (1._MK-2._MK*rho2+3./4.*rho2**2-1./12.*rho2**3+1./384.*rho2**4);
                factor = (1._MK-Qm*E)/r2;
                ! --------------------------------------------------------------------------------
                ! --- Mads Kernel (they are the same, it's just the fact that there is a 2 in the exponential that changes things)
                ! --------------------------------------------------------------------------------
            case(102) ! Second order Gaussian - Mads
                rho2   = r2/SmoothParam**2;
                E      = exp(-rho2/2._MK);  !Note the 2 in the Gaussian
                Qm     = 1._MK;
                factor = (1._MK-Qm*E)/r2;
            case(104) ! Fourth order Gaussian - Mads
                rho2   = r2/SmoothParam**2;
                E      = exp(-rho2/2._MK);  !Note the 2 in the Gaussian
                Qm     = 1._MK-1._MK/2._MK*rho2;
                factor = (1._MK-Qm*E)/r2;
            case(106)  ! Sixth Order - Gaussian - Mads
                rho2   = r2/SmoothParam**2;
                E      = exp(-rho2/2._MK);    !Note the 2 in the Gaussian
                Qm     = (1._MK-rho2+1._MK/8._MK*rho2**2);
                factor = (1._MK-Qm*E)/r2;
            case(108)  ! Eigth Order - Gaussian - Mads
                rho2   = r2/SmoothParam**2;
                E      = exp(-rho2/2._MK);                    !Note the 2 in the Gaussian
                Qm     = (1._MK-3._MK/2._MK*rho2+3._MK/8._MK*rho2**2-1._MK/48._MK*rho2**3);
                factor = (1._MK-Qm*E)/r2;
            case(110)  ! Tenth Order - Gaussian - Mads
                rho2   = r2/SmoothParam**2;
                E      = exp(-rho2/2._MK);                               !Note the 2 in the Gaussian
                Qm     = (1._MK-2._MK*rho2+3./4.*rho2**2-1./12.*rho2**3+1./384.*rho2**4);
                factor = (1._MK-Qm*E)/r2;
            case default 
                write(*,*)'2D vortex point: wrong smooth model'
                STOP -1
            end select 
            factor = factor * Intensity * twopi_inv;
            Uind(1)=-DeltaP(2)*factor
            Uind(2)= DeltaP(1)*factor
            Uind(3)=0.0_MK

        endif

    end subroutine

end module VortexPoint2D
