/*
* C function dlminf_nonplanarso.c -
*
* Calculates influence coefficient matrices for the nonplanar Doublet lattice
* Method.
*
* Inputs:
* xc: mc*nc matrix of x-coordinates of panel control points
* yc: mc*nc matrix of y-coordinates of panel control points
* zc: mc*nc matrix of z-coordinates of panel control points
* gammac: mc*nc matrix of dihedral angles of influenced panels
* xd: mv*(nv+1) matrix of x-coordinates of vertices of doublet lines
* yd: mv*(nv+1) matrix of y-coordinates of vertices of doublet lines
* zd: mv*(nv+1) matrix of z-coordinates of vertices of doublet lines
* xm: mv*nv matrix of x-coordinates of midpoints of doublet lines
* ym: mv*nv matrix of y-coordinates of midpoints of doublet lines
* zm: mv*nv matrix of z-coordinates of midpoints of doublet lines
* gammas: mv*nv matrix of dihedral angles of influencing panels
* Deltax: mv*nv matrix of mean chords of panels
* Deltay: (mv+1)*nv matrix of half-spans of panels
* params: 1*2 vector of parameters. The first element is the reduced
*         frequency k and the second the free stream Mach number.
* Areal: (mc*nc)*(mv*nv) matrix real part of influence coefficient matrix
* Aimag: (mc*nc)*(mv*nv) matrix imaginary part of influence coefficient
*        matrix
*
* Outputs: This code does not return any outputs. The results of its
*          calculations are stored in input matrices Areal and Aimag.
* NB: All matrices must have compatible dimensions; if they do not, C
*     will crash without warning.
*
 * Compile this code by typing at the terminal 
 * gcc -fPIC -shared -o dlminf_nonplanarso.so dlminf_nonplanarso.c
*
* This code is part of the SDPMflut Matlab distribution.
* Copyright (C) 2025 Grigorios Dimitriadis
*
* This program is free software: you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* This program is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
* Public License for more details.
*
* You should have received a copy of the GNU General Public License along
* with this program. If not, see <https://www.gnu.org/licenses/>.
*/

#include "math.h"
#include <stddef.h>

void I1I2(double IJ12ri[],double u1,double k1, double a[],double c)
{
    double coeff,denom,coeff1,coeff3,J1real,J1imag,J2real,J2imag;
    double I1real,I1imag,I2real,I2imag,sink1u1,cosk1u1;
    int n,np1;

    // Calculate sine and cosine of k1*u1
    sink1u1=sin(k1*u1);
    cosk1u1=cos(k1*u1);
    // Initialize J1 integral
    J1real=0.0;
    J1imag=0.0;
    // Initialize J2 integral
    J2real=0.0;
    J2imag=0.0;
    // Calculate J1 and J2 integrals
    for (n=0; n < 11; n++){
        np1=n+1;
        denom=(np1*np1*c*c+k1*k1);
        coeff=a[n]*exp(-np1*c*u1)/denom;
        J1real=J1real+coeff*np1*c;
        J1imag=J1imag-coeff*k1;
        coeff=coeff/denom;
        J2real=J2real+coeff*(np1*np1*c*c-k1*k1+np1*c*u1*denom);
        J2imag=J2imag-coeff*k1*(2*np1*c+u1*denom);
    }
    // Calculate I1 integral
    coeff=1.0-u1/sqrt(1.0+u1*u1)+k1*J1imag;
    I1real=cosk1u1*coeff-sink1u1*k1*J1real;
    I1imag= -cosk1u1*k1*J1real-sink1u1*coeff;
    // Calculate I2 integral
    coeff1=1-u1/sqrt(1+u1*u1);
    coeff3=u1/pow(1+u1*u1,1.5);
    coeff=(2.0*coeff1-coeff3+k1*J1imag+k1*k1*J2real);
    denom=k1*(u1*coeff1-J1real+k1*J2imag);
    I2real=1.0/3.0*(cosk1u1*coeff+sink1u1*denom);
    I2imag=1.0/3.0*(cosk1u1*denom-sink1u1*coeff);

    // Return integals in array IJ12ri
    IJ12ri[0]=I1real;
    IJ12ri[1]=I1imag;
    IJ12ri[2]=I2real;
    IJ12ri[3]=I2imag;
}

void dlminf_nonplanar(double *xc, double *yc, double *zc, double *xd, double *yd, double *zd, 
            double *xm, double *ym, double *zm, double *gammac, double *gammas,
			double *Deltax, double *Deltay, double *params, double *Areal, double *Aimag,
			int mc, int nc, int mv, int nv)
{
    double c,pi,a[11],beta2,beta,xdbar[2],ydbar[2],DeltaxJ,DeltayJ,LambdaJ;
    double k,Mach,xcIbar,ycIbar,zcIbar,gammaIbar,xbar[3],ybar[3],x0[3],y0[3],z0[3],zratio;
    double r[3],rbeta[3],u1[3],k1[3],I1real[3],I1imag[3],I2real[3],I2imag[3],IJ12ri[4],DeltayJdb2;
    double K1real,K1imag,P1real[3],P1imag[3],A1real,A1imag,B1real,B1imag,C1real,C1imag;
    double K2real,K2imag,P2real[3],P2imag[3],A2real,A2imag,B2real,B2imag,C2real,C2imag;
    double coeff,T1,T2star,num,denom,denom2,D1real,D1imag,D2real,D2imag,epsilon,aF,F,G,coeff1,coeff2;
    double I1minusreal,I1minusimag,I1zeroreal,I1zeroimag,I2minusreal,I2minusimag,I2zeroreal,I2zeroimag;
    double K,L,M;
    int    mcnc,mvnv,nvp1,ic,jc,iv,jv,ik;

    /* Determine sizes of input arrays */
    mcnc=mc*nc;
    mvnv=mv*nv;
    nvp1=nv+1;
    
    /* Assign value of Mach */
    k= *params;
    Mach= *(params+1);
    pi=3.14159265358979;  /* value of pi */

    beta2=1-Mach*Mach;
    beta=sqrt(beta2);
    /* Laschka approximation coefficients */
    c=0.372;
    a[0]= 0.24186198;
    a[1]= -2.7918027;
    a[2]= 24.991079;
    a[3]= -111.59196;
    a[4]= 271.43549;
    a[5]= -305.75288;
    a[6]= -41.183630;
    a[7]= 545.98537;
    a[8]= -644.78155;
    a[9]= 328.72755;
    a[10]= -64.279511;

    /* Create influence coefficient matrix */
    /* Cycle through influenced control points */
    for (ic=0; ic < mc; ic++) {
        for (jc=0; jc < nc; jc++) {
            /* Cycle through influencing doublet panels */
            for (iv=0; iv < mv; iv++) {
                for (jv=0; jv < nv; jv++) {
                    // Calculate endpoints of Jth doublet line in Jth panel coordinates
                    xdbar[0]= *(xd+iv*nvp1+jv)- *(xm+iv*nv+jv);
                    xdbar[1]= *(xd+iv*nvp1+jv+1)- *(xm+iv*nv+jv);
                    ydbar[0]= (*(yd+iv*nvp1+jv)- *(ym+iv*nv+jv))*cos(*(gammas+iv*nv+jv))+(*(zd+iv*nvp1+jv)- *(zm+iv*nv+jv))*sin(*(gammas+iv*nv+jv));
                    ydbar[1]= (*(yd+iv*nvp1+jv+1)- *(ym+iv*nv+jv))*cos(*(gammas+iv*nv+jv))+(*(zd+iv*nvp1+jv+1)- *(zm+iv*nv+jv))*sin(*(gammas+iv*nv+jv));
                    // Read mean chordwise length of Jth panel
                    DeltaxJ= *(Deltax+iv*nv+jv);
                    // Read mean spanwise length of Jth panel
                    DeltayJ=2* *(Deltay+iv*nv+jv);
                    // Half mean spanwise length of Jth panel
                    DeltayJdb2=DeltayJ/2.0;
                    // Calculate sweep of Jth doublet line
                    LambdaJ=atan2(xdbar[1]-xdbar[0],ydbar[1]-ydbar[0]);
                    // Calculate coordinates of Ith control point in Jth panel coordinates
                    xcIbar= *(xc+ic*nc+jc)- *(xm+iv*nv+jv);
                    ycIbar=(*(yc+ic*nc+jc)- *(ym+iv*nv+jv))*cos(*(gammas+iv*nv+jv))+(*(zc+ic*nc+jc)- *(zm+iv*nv+jv))*sin(*(gammas+iv*nv+jv));
                    zcIbar=(*(zc+ic*nc+jc)- *(zm+iv*nv+jv))*cos(*(gammas+iv*nv+jv))-(*(yc+ic*nc+jc)- *(ym+iv*nv+jv))*sin(*(gammas+iv*nv+jv));
                    // Calculate dihedral angle of Ith control point in Jth panel coordinates
                    gammaIbar= *(gammac+ic*nc+jc)- *(gammas+iv*nv+jv);
                    // Calculate spanwise coordinates of Jth doublet line: two ends and midpoint
                    ybar[0]= -DeltayJdb2;
                    ybar[1]= 0.0;
                    ybar[2]= DeltayJdb2;
                   
                    for (ik=0; ik<3; ik++){
                        // Calculate chordwise coordinates of Jth doublet line
                        xbar[ik]=ybar[ik]*tan(LambdaJ);
                        // Calculate distance between Ith control point and points on doublet line
                        x0[ik]=xcIbar-xbar[ik];
                        y0[ik]=ycIbar-ybar[ik];
                        z0[ik]=zcIbar;
                        // Calculate r distance
                        r[ik]=sqrt(y0[ik]*y0[ik]+z0[ik]*z0[ik]);
                        // Remove zeros from r
                        if (r[ik] < 1e-12){
                            r[ik]=1e-12;
                        }
                        // Calculate rbeta distance
                        rbeta[ik]=sqrt(x0[ik]*x0[ik]+beta2*r[ik]*r[ik]);
                        // Calculate u1
                        u1[ik]=(-x0[ik]+Mach*rbeta[ik])/r[ik]/beta2;
                        // Calculate k1
                        k1[ik]=k*r[ik];
                    }

                    // Calculate integral of Kernel function
                    for (ik=0; ik<3; ik++) {
                        if (u1[ik] >= 0.0){
                            I1I2(IJ12ri,u1[ik],k1[ik],a,c);
                            I1real[ik]=IJ12ri[0];
                            I1imag[ik]=IJ12ri[1];
                            I2real[ik]=IJ12ri[2];
                            I2imag[ik]=IJ12ri[3];
                        }else{
                            I1I2(IJ12ri,-u1[ik],k1[ik],a,c);
                            I1minusreal=IJ12ri[0];
                            I1minusimag=IJ12ri[1];
                            I2minusreal=IJ12ri[2];
                            I2minusimag=IJ12ri[3];
                            I1I2(IJ12ri,0.0,k1[ik],a,c);
                            I1zeroreal=IJ12ri[0];
                            I1zeroimag=IJ12ri[1];
                            I2zeroreal=IJ12ri[2];
                            I2zeroimag=IJ12ri[3];
                            I1real[ik]=2.0*I1zeroreal-I1minusreal;
                            I1imag[ik]=I1minusimag;
                            I2real[ik]=2.0*I2zeroreal-I2minusreal;
                            I2imag[ik]=I2minusimag;
                        }
                    }

                    // Calculate T1, K1, P1 and numerator of surface integral
                    T1=cos(gammaIbar);
                    for (ik=0; ik<3; ik++){
                        coeff=Mach*r[ik]/rbeta[ik]/sqrt(1+u1[ik]*u1[ik]);
                        K1real=I1real[ik]+coeff*cos(k1[ik]*u1[ik]);
                        K1imag=I1imag[ik]-coeff*sin(k1[ik]*u1[ik]);
                        P1real[ik]=(cos(k*x0[ik])*K1real+sin(k*x0[ik])*K1imag)*T1;
                        P1imag[ik]=(cos(k*x0[ik])*K1imag-sin(k*x0[ik])*K1real)*T1;
                    }

                    // Calculate coefficients of quadratic approximation of numerator
                    A1real=(P1real[0]-2*P1real[1]+P1real[2])/2.0/DeltayJdb2/DeltayJdb2;
                    A1imag=(P1imag[0]-2*P1imag[1]+P1imag[2])/2.0/DeltayJdb2/DeltayJdb2;
                    B1real=(P1real[2]-P1real[0])/DeltayJ;
                    B1imag=(P1imag[2]-P1imag[0])/DeltayJ;
                    C1real=P1real[1];
                    C1imag=P1imag[1];

                    // Calculate logarithmic term
                    num=(ycIbar-DeltayJdb2)*(ycIbar-DeltayJdb2)+zcIbar*zcIbar;
                    denom=(ycIbar+DeltayJdb2)*(ycIbar+DeltayJdb2)+zcIbar*zcIbar;
                    G=log(num/denom);

                    // Calculate planarity ratio
                    zratio=fabs(zcIbar)/DeltayJdb2;

                    if (zratio < 0.001){
                        // Planar case
                        // Calculate F
                        F=DeltayJ/(ycIbar*ycIbar-DeltayJdb2*DeltayJdb2);
                        //  Calculate planar upwash factor
                        D1real=DeltaxJ*((ycIbar*ycIbar*A1real+ycIbar*B1real+C1real)*F \
                                        +(B1real/2.0+ycIbar*A1real)*G+DeltayJ*A1real)/8.0/pi;
                        D1imag=DeltaxJ*((ycIbar*ycIbar*A1imag+ycIbar*B1imag+C1imag)*F \
                                        +(B1imag/2.0+ycIbar*A1imag)*G+DeltayJ*A1imag)/8.0/pi;
                        // Set nonplanar upwash factor to zero
                        D2real=0.0;
                        D2imag=0.0;
                    }else{
                        // Nonplanar case
                        // Calculate term to decide if we will use an approximation for F
                        denom=ycIbar*ycIbar+zcIbar*zcIbar-DeltayJdb2*DeltayJdb2;
                        epsilon=DeltayJ*fabs(zcIbar)/denom;
                        if (epsilon <= 0.3){
                            aF=0;
                            for (ik=2; ik<8; ik++) {
                                aF+= pow(-1.0,ik)/(2*ik-1)*pow(epsilon,2*ik-4);
                            }
                            aF=pow(DeltayJ,4.0)/4.0/denom/denom*aF;
                            F=DeltayJ/denom*(1-aF*zcIbar*zcIbar/(DeltayJdb2*DeltayJdb2));
                        }else{
                            F=1/fabs(zcIbar)*atan2(DeltayJ*fabs(zcIbar),denom);
                        }
                        // Calculate planar upwash factor
                        D1real=DeltaxJ*(((ycIbar*ycIbar-zcIbar*zcIbar)*A1real+ycIbar*B1real+C1real)*F \
                                        +(B1real/2+ycIbar*A1real)*G+DeltayJ*A1real)/8.0/pi;
                        D1imag=DeltaxJ*(((ycIbar*ycIbar-zcIbar*zcIbar)*A1imag+ycIbar*B1imag+C1imag)*F \
                                        +(B1imag/2+ycIbar*A1imag)*G+DeltayJ*A1imag)/8.0/pi;

                        // Calculate T2star, K2, P2 and numerator of surface integral
                        for (ik=0; ik<3; ik++){
                            T2star=(z0[ik]*cos(gammaIbar)-y0[ik]*sin(gammaIbar))*z0[ik];
                            coeff=Mach*r[ik]/rbeta[ik]*((1+u1[ik]*u1[ik])*beta2*r[ik]*r[ik]/rbeta[ik]/rbeta[ik] \
                                                        +2.0+Mach*r[ik]*u1[ik]/rbeta[ik])/pow(1+u1[ik]*u1[ik],1.5);
                            num=k1[ik]*Mach*Mach*r[ik]*r[ik]/rbeta[ik]/rbeta[ik]/sqrt(1+u1[ik]*u1[ik]);
                            K2real= -3.0*I2real[ik]-num*sin(k1[ik]*u1[ik])-coeff*cos(k1[ik]*u1[ik]);
                            K2imag= -3.0*I2imag[ik]-num*cos(k1[ik]*u1[ik])+coeff*sin(k1[ik]*u1[ik]);
                            P2real[ik]=(cos(k*x0[ik])*K2real+sin(k*x0[ik])*K2imag)*T2star;
                            P2imag[ik]=(cos(k*x0[ik])*K2imag-sin(k*x0[ik])*K2real)*T2star;
                        }

                        // Calculate coefficients of quadratic approximation of numerator
                        A2real=(P2real[0]-2*P2real[1]+P2real[2])/2.0/DeltayJdb2/DeltayJdb2;
                        A2imag=(P2imag[0]-2*P2imag[1]+P2imag[2])/2.0/DeltayJdb2/DeltayJdb2;
                        B2real=(P2real[2]-P2real[0])/DeltayJ;
                        B2imag=(P2imag[2]-P2imag[0])/DeltayJ;
                        C2real=P2real[1];
                        C2imag=P2imag[1];

                        denom=((ycIbar+DeltayJdb2)*(ycIbar+DeltayJdb2)+zcIbar*zcIbar);
                        denom2=((ycIbar-DeltayJdb2)*(ycIbar-DeltayJdb2)+zcIbar*zcIbar);
                        if (1.0/fabs(epsilon) <= 0.1){
                            // Calculate K, L, M terms
                            coeff=ycIbar*ycIbar+zcIbar*zcIbar;
                            num=ycIbar*ycIbar-zcIbar*zcIbar;
                            K=coeff*F+(coeff*ycIbar+num*DeltayJdb2)/denom-(coeff*ycIbar-num*DeltayJdb2)/denom2;
                            L=ycIbar*F+(coeff+ycIbar*DeltayJdb2)/denom-(coeff-ycIbar*DeltayJdb2)/denom2;
                            M=F+(ycIbar+DeltayJdb2)/denom-(ycIbar-DeltayJdb2)/denom2;
                            // Calculate nonplanar upwash factor
                            D2real=DeltaxJ/2.0/zcIbar/zcIbar*(K*A2real+L*B2real+M*C2real)/8.0/pi;
                            D2imag=DeltaxJ/2.0/zcIbar/zcIbar*(K*A2imag+L*B2imag+M*C2imag)/8.0/pi;
                        }else{
                            // Calculate aF
                            coeff=ycIbar*ycIbar+zcIbar*zcIbar-DeltayJdb2*DeltayJdb2;
                            num=ycIbar*ycIbar+zcIbar*zcIbar+DeltayJdb2*DeltayJdb2;
                            aF=DeltayJdb2*DeltayJdb2/zcIbar/zcIbar*(1-coeff*F/DeltayJ);
                            coeff1=DeltayJdb2*DeltaxJ/coeff;
                            coeff2=1/denom/denom2;
                            K=coeff2*2.0*num*DeltayJdb2*DeltayJdb2-aF/DeltayJdb2/DeltayJdb2*(ycIbar*ycIbar+zcIbar*zcIbar);
                            L=coeff2*ycIbar*DeltayJ*DeltayJ-aF/DeltayJdb2/DeltayJdb2*ycIbar;
                            M=coeff2*2.0*num-aF/DeltayJdb2/DeltayJdb2;
                            D2real=coeff1*(K*A2real+L*B2real+M*C2real)/8.0/pi;
                            D2imag=coeff1*(K*A2imag+L*B2imag+M*C2imag)/8.0/pi;
                        }
                    }
                    // Total influence coefficient matrix
                    *(Areal+(ic*nc+jc)*mvnv+(iv*nv+jv))=D1real+D2real;
                    *(Aimag+(ic*nc+jc)*mvnv+(iv*nv+jv))=D1imag+D2imag;
                }
            }
        }
    }
}
