Module invariants_mod
    implicit none 
    contains
    !----------------------------------------------------
    subroutine invariants(a,xs,ys,xb,yb,b,inv)
        implicit none 
        !include "nnparams.incl"
        double precision, intent(in) :: a, xs, ys, xb, yb, b
        double precision, intent(out) :: inv(4)
        double precision:: invar(24)
        complex(8) :: q1, q2
        LOGICAL,PARAMETER:: debg =.false.
        integer :: i 
        ! express the coordinate in complex

        q1 = dcmplx(xs, ys)
        q2 = dcmplx(xb, yb)

        ! compute the invariants
        invar(24) = a 
        invar(23) =b**2 

        ! INVARIANTS OF KIND II
        !------------------------

        invar(1) = dreal( q1 * conjg(q1) ) ! r11
        invar(2) = dreal( q1 * conjg(q2) ) ! r12
        invar(3) = dreal( q2 * conjg(q2) ) ! r22
        invar(4) = (dimag(q1 * conjg(q2)) )**2 ! rho 12**2  


        !INVATIANTS OF KIND III 
        !------------------------

        invar(5) = dreal( q1 * q1 * q1 ) ! r111
        invar(6) = dreal( q1 * q1 * q2 ) ! r112
        invar(7) = dreal( q1 * q2 * q2 ) ! r122
        invar(8) = dreal( q2 * q2 * q2 ) ! r222
        invar(9) = (dimag( q1 * q1 * q1 ))**2 ! rho111**2
        invar(10) = (dimag( q1 * q1 * q2 ))**2 ! rho112 **2
        invar(11) = (dimag( q1 * q2 * q2 ))**2 ! rho122**2 
        invar(12) = (dimag( q2 * q2 * q2 ))**2 ! rho222

        ! INVARIANTS OF KIND IV 
        !-------------------------

        invar(13) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q1 ))
        invar(14) = (dimag( q1 * conjg(q2)) * dimag( q1 * q1 * q2 ))
        invar(15) = (dimag( q1 * conjg(q2)) * dimag( q1 * q2 * q2 ))
        invar(16) = (dimag( q1 * conjg(q2)) * dimag( q2 * q2 * q2 ))

        ! INVARIANTS OF KIND V
        !----------------------

        invar(17) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q1 * q2 ))
        invar(18) = (dimag( q1 * q1 * q1 ) * dimag( q1 * q2 * q2 ))
        invar(19) = (dimag( q1 * q1 * q1 ) * dimag( q2 * q2 * q2 ))
        invar(20) = (dimag( q1 * q1 * q2 ) * dimag( q1 * q2 * q2 ))
        invar(21) = (dimag( q1 * q1 * q2 ) * dimag( q2 * q2 * q2 ))
        invar(22) = (dimag( q1 * q2 * q2 ) * dimag( q2 * q2 * q2 ))

        ! the only non zero invariant for bend pure cuts

        inv(1) = invar(1)
        inv(2) = invar(5)
        inv(3) = invar(9)
        inv(4) = invar(23)
        if (debg) then 
            write(14,"(A,*(f10.5))")"Invar II", (invar(i),i=1,4)
            write(14,"(A,*(f10.5))") "Invar III", (invar(i),i=5,12)
            write(14,"(A,*(f10.5))")"Invar IV", (invar(i),i=13,16)
            write(14,"(A,*(f10.5))")"Invar V", (invar(i),i=17,22)

            write(14,*)"THE INPUT COORDINATE IN COMPLEX REPRES"
            write(14,*)"---------------------------------------"
            write(14,*)"xs =",dreal(q1), "ys=",dimag(q1)
        endif     
    ! modify the invariants to only consider few of them 
    ! 
    !invar(13:22)=0.0d0
    end subroutine invariants
end module invariants_mod