using SpecialFunctions; # need to evaluate the Hankel function (Green's function)
using Luxor;  # this package is needed to indentify the points inside/outside the curve
using Plots;
#--------------------------------------------------------------------------------------- 
struct Curve
    name::String
    N::Int64               # number of nodes
    x::Array{Float64,2}    # nodes on the curve generated by a parametrization
    dx::Array{Float64,2}   # 1st derivative of the curve parametrization at the nodes 
    d2x::Array{Float64,2}  # 2nd derivative of the curve parametrization at the nodes 
    d3x::Array{Float64,2}  # 3rd derivative of the curve parametrization at the nodes     
    τ::Array{Float64,1}    # norm of the derivative of the curve parametrization at the nodes
    ν::Array{Float64,2}    # unit normal at the nodes
    Δs::Float64            # discretization size in parameter space (~2*pi/N)
    s::Array{Float64,1}    # discrete values of the parameter 
    ds::Array{Float64,1}   # derivative of the sigmoid change of variable (for curves with corners)
    R::Array{Float64,2}    # Nystrom weights    
    param::Function        # Curve parametrization
    PD:: Int64             # Polynomial degree in the mesh refinement
end
#---------------------------------------------------------------------------------------
function SkeletonCctor(ver,dir,M::Array{Int64,1},p=4)
    # ver: coordinates of the vertices
    # dir: orientation of the segment: 
    #       dir[1]:index of starting vertex
    #       dir[2]:index of ending vertex
    # M: roughly half of the number of nodes in the segment
    # p: polinomial degree of the sigmoid change of variables

    Np    = size(dir,1);
    parts = Vector{Curve}(undef,Np)

    Mvec = Array{Int64,1}(undef,Np)
    if length(M)==1
        Mvec .= M
    else
        Mvec = M
    end
    
    Npts = 0
    for n=1:Np
        st = dir[n,1];ed = dir[n,2];
        parts[n] = lineSegment(Mvec[n],ver[st,:],ver[ed,:],p);                
        Npts = Npts+parts[n].N
    end

    plt = plot(parts[1].x[:,1],parts[1].x[:,2]);
    for n =2:Np
        plt = plot!(parts[n].x[:,1],parts[n].x[:,2]);
    end
    plt = plot!(plt,legend=false,aspect_ratio=:equal,frame=:box)
    # return
    return s = (parts=parts,vertices=ver,orientation=dir,Nparts=Np,plt=plt,Npts = Npts)

end
#--------------------------------------------------------------------------------------- 
function NystromWeights(M)
    R = zeros(2*M,2*M);
        for p=1:2*M            
            tp = pi/M*(p-1);
            for j=1:2*M
                tj = pi/M*(j-1);
                R[p,j]=-2*pi/M*sum((ones(M-1)./collect(1:M-1)).*cos.(collect(1:M-1)*(tp-tj)))-pi/(M^2)*cos(M*(tp-tj));               
            end            
        end
    return R
end

#--------------------------------------------------------------------------------------- 

function lineSegment(M::Int64,xs,xe,p=4)
    Δs = pi/M
    t = Δs*collect(1:2*M-1)
    N = 2*M-1
    (s,ds) = polySigmoid(t,p)

    param(s) = ([xs[1].+(xe[1]-xs[1])/2/pi*s xs[2].+(xe[2]-xs[2])/2/pi*s],1/2/pi*[(xe[1]-xs[1])*size(s,1) (xe[2]-xs[2])*size(s,1)])
    
    # Segment
    x   = [xs[1].+(xe[1]-xs[1])/2/pi*s xs[2].+(xe[2]-xs[2])/2/pi*s]
    dx  = 1/2/pi*[(xe[1]-xs[1])*ones(N,1) (xe[2]-xs[2])*ones(N,1)]
    d2x = zeros(N,2)
    d3x = zeros(N,2)
    
    τ  = sqrt.(dx[:,1].^2+dx[:,2].^2);
    ν = [dx[:,2]./τ -dx[:,1]./τ] 

    c = Curve("line",N,x,dx,d2x,d3x,τ,ν,Δs,s,ds,NystromWeights(M),param,p)

    return c

end

# --------------------------------------------------------------------------------------- 

function kite(M::Integer,rad=1,loc=[0 0])
    Δs = pi/M
    s = Δs*collect(0:2*M-1);
    N = 2*M
    param(s) =(loc.+rad*[cos.(s)+0.65*cos.(2*s).-0.65 1.5*sin.(s)],rad*[-sin.(s)-1.3*sin.(2*s) 1.5*cos.(s)])
    # kite    
    x  =loc.+rad*[cos.(s)+0.65*cos.(2*s).-0.65 1.5*sin.(s)];
    dx = rad*[-sin.(s)-1.3*sin.(2*s) 1.5*cos.(s)];
    d2x = rad*[-cos.(s)-2.6*cos.(2*s) -1.5*sin.(s)];
    d3x = rad*[sin.(s)+5.2*sin.(2*s) -1.5*cos.(s)];

    τ  = sqrt.(dx[:,1].^2+dx[:,2].^2);
    ν = [dx[:,2]./τ -dx[:,1]./τ] 
    ds = ones(N);

    c = Curve("kite",N,x,dx,d2x,d3x,τ,ν,Δs,s,ds,NystromWeights(M),param,0);

    return c    
end
# ---------------------------------------------------------------------------------------
function matricesSkeleton(prtId::Int,sk,kVec,coeff=[1.0 -1;1 -1;1 -1;1 -1])
    
    prt = sk.parts[prtId]

    if length(kVec)==1

        k=kVec

        if 1==prtId
            matAux = matricesNear(k,prt)                               
        else
            matAux = matricesFar(k,prt,sk.parts[1])
        end

        S  = matAux.S;
        K  = matAux.K;
        Kp = matAux.Kp;
        T  = matAux.T;

        Np = sk.Nparts;

        for n=2:Np
            if n==prtId
                matAux = matricesNear(k,prt)                            
            else
                matAux = matricesFar(k,prt,sk.parts[n])
            end
            #println(size(matAux.S))

            S  =[S matAux.S];
            K  =[K matAux.K];
            Kp =[Kp matAux.Kp];
            T  =[T matAux.T];

        end

        return mat = (S=S,K=K,Kp=Kp,T=T)

    elseif length(kVec)==2
        
        matAux1 = matricesSkeleton(prtId,sk,kVec[1])
        matAux2 = matricesSkeleton(prtId,sk,kVec[2])

        S  = coeff[1,1]*matAux1.S  + coeff[1,2]*matAux2.S
        K  = coeff[2,1]*matAux1.K  + coeff[2,2]*matAux2.K
        Kp = coeff[3,1]*matAux1.Kp + coeff[3,2]*matAux2.Kp
        T  = coeff[4,1]*matAux1.T  + coeff[4,2]*matAux2.T

        return mat = (S=S,K=K,Kp=Kp,T=T)

    end
    
end
# ---------------------------------------------------------------------------------------
function matricesNear(k,prt)
        
    τ = prt.τ
    ν = prt.ν
    s = prt.s
    ds= prt.ds
    R = prt.R

    # R = NystromWeights(prt.N)

    N = prt.N
    Δs = prt.Δs

    x1 = prt.x[:,1];     x2 = prt.x[:,2]
    dx1 = prt.dx[:,1];   dx2 = prt.dx[:,2]
    d2x1 = prt.d2x[:,1]; d2x2 = prt.d2x[:,2]
    d3x1 = prt.d3x[:,1]; d3x2 = prt.d3x[:,2]
    n1 = ν[:,1];n2=ν[:,2]
    
    
    S  = zeros(N,N)+im*zeros(N,N)
    K  = zeros(N,N)+im*zeros(N,N)
    Kp = zeros(N,N)+im*zeros(N,N)
    T  = zeros(N,N)+im*zeros(N,N)
    
    psi = -0.577215664901532;
    
    ii = floor(N/2)
    for i=1:N
        for j=1:N

            r = sqrt((x1[i]-x1[j])^2+(x2[i]-x2[j])^2);            

            J0 = besselj0(k*r); J1 = besselj1(k*r)
            H0 = hankelh1(0,k*r); H1 = hankelh1(1,k*r)
            
            lSin = log(4*(sin(Δs*(i-j)/2))^2)

            if i!=j

                S1 = (-1/(4*pi))*J0*τ[j];                
                S2 = im/4*H0*τ[j]-S1*lSin

                K1 = -k/(4*pi)*J1/r*(dx2[j]*(x1[i]-x1[j])-dx1[j]*(x2[i]-x2[j]));        
                K2 = im*k/4*H1/r*((x1[i]-x1[j])*dx2[j]-(x2[i]-x2[j])*dx1[j])-K1*lSin

                Kp1 = k/(4*pi)*J1/r*(dx2[i]*(x1[i]-x1[j])-dx1[i]*(x2[i]-x2[j]))*τ[j]/τ[i]
                Kp2 = im*k/4*H1/r*((x1[j]-x1[i])*dx2[i]-(x2[j]-x2[i])*dx1[i])*τ[j]/τ[i]-Kp1*lSin
                
                Nt = ((x1[i]-x1[j])*dx1[j]+(x2[i]-x2[j])*dx2[j])*((x1[i]-x1[j])*dx1[i]+(x2[i]-x2[j])*dx2[i])/r^2

                NS1 = -1/2/pi*Nt*(k^2*J0-2*k*J1/r)-k/2/pi*(dx1[i]*dx1[j]+dx2[i]*dx2[j])/r*J1

                NS = im/2*Nt*(k^2*H0-2*k*H1/r)+im*k/2*(dx1[i]*dx1[j]+dx2[i]*dx2[j])/r*H1+1/4/pi/(sin((s[i]-s[j])/2))^2

                NS2 = NS-NS1*lSin

            else

                S1 = (-1/(4*pi))*J0*τ[j];                            
                S2  =(im/4+psi/2/pi-1/4/pi*log(k^2/4*τ[j]^2))*τ[j]+2*log(ds[j])*S1

                K1 = 0.0
                K2 = -1/(4*pi)*(dx1[j]*d2x2[j]-dx2[j]*d2x1[j])/τ[j]^2+2*log(ds[j])*K1
                
                Kp1 = K1; Kp2 = K2

                NS1 = -k^2/4/pi*τ[i]^2

                NS2 = k^2/4/pi*(pi*im-1+2*psi-2*log(k*τ[i]/2))*τ[i]^2+1/12/pi+            
                    +1/2/pi*(dx1[i]*d2x1[i]+dx2[i]*d2x2[i])^2/τ[i]^4 +
                    -1/4/pi*(d2x1[i]^2+d2x2[i]^2)/τ[i]^2 +
                    -1/6/pi*(dx1[i]*d3x1[i]+dx2[i]*d3x2[i])/τ[i]^2

                NS2 = NS2 + 2*log(ds[j])*NS1

            end 

            S[i,j] = (R[i,j]*S1 + Δs*S2)*ds[j];
            K[i,j] = (R[i,j]*K1 + Δs*K2)*ds[j];
            Kp[i,j] = (R[i,j]*Kp1 + Δs*Kp2)*ds[j];
            
            T[i,j] = -(R[i,j]*NS1 + Δs*NS2)*ds[j]/2/τ[i];
            T[i,j] = k^2*S[i,j]*(ν[i,1]*ν[j,1]+ν[i,2]*ν[j,2]) + T[i,j];
           
        end

    end
    

    mat = (S=S,K=K,Kp=Kp,T=T);

    plot(plt)
    return mat

end
# --------------------------------------------------------------------------------------- 
function matricesFar(k,prtEvl,prtInt)
    
    τI = prtInt.τ
    τE = prtEvl.τ

    νI = prtInt.ν
    νE = prtEvl.ν
    
    ds = prtInt.ds
    Δs = prtInt.Δs

    NI = prtInt.N
    NE = prtEvl.N

    x1I = prtInt.x[:,1]
    x2I = prtInt.x[:,2]

    x1E = prtEvl.x[:,1]
    x2E = prtEvl.x[:,2]

    dx1I = prtInt.dx[:,1]
    dx2I = prtInt.dx[:,2]

    dx1E = prtEvl.dx[:,1]
    dx2E = prtEvl.dx[:,2]

    S = zeros(NE,NI)+im*zeros(NE,NI)
    K = zeros(NE,NI)+im*zeros(NE,NI)
    Kp = zeros(NE,NI)+im*zeros(NE,NI)
    T = zeros(NE,NI)+im*zeros(NE,NI)
        
    for i=1:NE,j=1:NI

        r = sqrt((x1E[i]-x1I[j])^2+(x2E[i]-x2I[j])^2)

        H0 = hankelh1(0,k*r)
        H1 = hankelh1(1,k*r)
        
        S2 = im/4*H0*τI[j];
        K2 = im*k/4*H1/r*((x1E[i]-x1I[j])*νI[j,1]+(x2E[i]-x2I[j])*νI[j,2])*τI[j]

        Kp2 = im*k/4*H1/r*((x1I[j]-x1E[i])*νE[i,1]+(x2I[j]-x2E[i])*νE[i,2])*τI[j]

        # single layer 
        S[i,j] = Δs*S2*ds[j];

        # double layer
        K[i,j] = Δs*K2*ds[j];        

        # adjoint double layer operator
        Kp[i,j] = Δs*Kp2*ds[j];

        # hipersingular
        # A11 = (im/4)*k^2*(x1E[i]-x1I[j])^2*H0/r^2 +
        #        (im/4)*k*(r^2-2*(x1E[i]-x1I[j])^2)*H1/r^3

        # A12 = -(im/2)*k*(x1E[i]-x1I[j])*(x2E[i]-x2I[j])*H1/r^3 +
        #        (1i/4)*k^2*(x1E[i]-x1I[j])*(x2E[i]-x2I[j])*H0/r^2

        # A21 = -(1i/2)*k*(x1E[i]-x1I[j]).*(x2E[i]-x2I[j])*H1/r^3 +
        #        (1i/4)*k^2*(x1E[i]-x1I[j])*(x2E[i]-x2I[j])*H0/r^2

        # A22 = (im/4)*k^2*(x2E[i]-x2I[j])^2*H0/r^2 +
        #       (im/4)*k*(r^2-2*(x2E[i]-x2I[j])^2)*H1/r^3
        
        # T[i,j] = (νE[i,1]*(A11*νI[j,1]+ A12*νI[j,2]) + νE[i,2]*(A21*νI[j,1] + A22*νI[j,2]))*Δs*τI[j]*ds[j]

        R1 = x1E[i]-x1I[j];R2 = x2E[i]-x2I[j];

        A11 = (im/4)*k^2*R1^2*H0/r^2+(im/4)*k*(r^2-2*R1^2)*H1/r^3

        A12 = -(im/2)*k*R1*R2*H1/r^3+(im/4)*k^2*R1*R2*H0/r^2

        A21 = -(im/2)*k*R1*R2*H1/r^3+(im/4)*k^2*R1*R2*H0/r^2

        A22 = (im/4)*k^2*R2^2*H0/r^2+(im/4)*k*(r^2-2*R2^2)*H1/r^3
    
        n1I =νI[j,1];n2I =νI[j,2]; 
        n1E =νE[i,1];n2E =νE[i,2]; 

        T[i,j] = (n1E*(A11*n1I + A12*n2I) + n2E*(A21*n1I + A22*n2I))*Δs*ds[j]*τI[j]
    

    end

    mat = (S=S,K=K,Kp=Kp,T=T);

    return mat

end
# --------------------------------------------------------------------------------------- 
function potentials(k,x,prt)
    
    Npe = size(x,1)
    Npi = prt.N

    SL = Array{Complex{Float64},2}(undef,Npe,Npi)
    DL = Array{Complex{Float64},2}(undef,Npe,Npi)
                
    x1 = x[:,1];x2 = x[:,2];
            
    xp1 = prt.x[:,1];xp2 = prt.x[:,2];
        
    # dxp1 = prt.dx[:,1];dxp2 = prt.dx[:,2];

    ν1 = prt.ν[:,1];ν2 = prt.ν[:,2]

    ds = prt.ds;Δs = prt.Δs
        
    tol  = 1e-4;

    for i = 1:Npe
        for j=1:Npi
            r = sqrt((x1[i]-xp1[j])^2+(x2[i]-xp2[j])^2)
            if r<tol
                r=tol
            end
            SL[i,j] = im/4*hankelh1(0,k*r)*prt.τ[j]*prt.ds[j]*Δs
            DL[i,j] = im*k/4*hankelh1(1,k*r)/r*((x1[i]-xp1[j])*ν1[j]+(x2[i]-xp2[j])*ν2[j])*prt.τ[j]*prt.ds[j]*Δs
        end
    end

    pot = (SL=SL,DL=DL);
    return pot
end

#--------------------------------------------------------
function potEval(φ,k,lims,h,crv,coeff)
        
    xLim = lims[1:2]
    yLim = lims[3:4]

    xa = xLim[1]:h:xLim[2]
    ya = yLim[1]:h:yLim[2]

    Nx = size(xa,1)
    Ny = size(ya,1)

    X = repeat(reshape(xa, 1, :), Ny, 1)
    Y = repeat(ya, 1, Nx)

    pts = [reshape(X,Nx*Ny,1) reshape(Y,Nx*Ny,1)];

    pot = potentials(k,pts,crv);

    Pφ = coeff[1]*pot.DL*φ+coeff[2]*pot.SL*φ;
    
    Pφ_in  = im*zeros(Nx*Ny,1)
    Pφ_out = im*zeros(Nx*Ny,1)
    
    # Find points inside the curve crv
    polygon = SVector.(crv.x[:,1],crv.x[:,2])
    points = vec(SVector.(xa',ya))
    inside = [isinside(p, polygon; allowonedge=true) for p in points]
    
    for j=1:Nx*Ny
        if inside[j]==true
            Pφ_in[j] = Pφ[j]
            Pφ_out[j] = NaN            
        else
            Pφ_out[j] = Pφ[j]
            Pφ_in[j] = NaN            
        end
    end
    
    potPlot=(In=reshape(Pφ_in,Ny,Nx),Out=reshape(Pφ_out,Ny,Nx),x=xa,y=ya,X=X,Y=Y,Nx=Nx,Ny=Ny);

    return potPlot
end

#--------------------------------------------------------
function solPlot(sol,crv,opt,reg,cscale=(-1,1));
    if opt=="real"
        if reg=="interior"
            plt = heatmap(sol.x, sol.y,real.(sol.In),color=:RdBu,clim=cscale)
        elseif reg=="exterior"
            plt = heatmap(sol.x, sol.y,real.(sol.Out),color=:RdBu,clim=cscale)
        end
    elseif opt=="abs"
        if reg=="interior"
            plt = heatmap(sol.x, sol.y,abs.(sol.In),clim=cscale)
        elseif reg=="exterior" 
            plt = heatmap(sol.x, sol.y,abs.(sol.Out),clim=cscale)
        end
    end

    pX = [crv.x[:,1];crv.x[1,1]];pY = [crv.x[:,2];crv.x[1,2]]  
    plt = plot!(pX,pY,aspect_ratio=:equal,lw=1.5,linecolor=:black,legend=false,frame=:box)
    
    return plt

end

#--------------------------------------------------------
function polySigmoid(s,p);
    #  Change of variable that cancels the singularity at the end points
    v(s,p) = (1/p-1/2)*((pi.-s)/pi).^3 .+1/p*(s.-pi)/pi .+ 1/2;
    dv(s,p)= -(3*(1/p-1/2))*(pi.-s).^2/pi^3 .+1/(p*pi)

    v0 = v(s,p)
    v1 = v(2*pi.-s,p)

    w = 2*pi*v0.^p./(v0.^p.+v1.^p)
    
    dv0 = dv(s,p);
    dv1 = dv(2*pi.-s,p);
    
    dw = 2*p*pi*(v0.^(p-1).*dv0.*v1.^p+v0.^p.*v1.^(p-1).*dv1)./(v0.^p+v1.^p).^2;

    w0  = (w=w,dw=dw)
    return w0;
end
#--------------------------------------------------------    
