/* Helium Peng-Robinson and compressible-orifice kernels.
 * Ported from the project's Python physical equations; validated independently
 * against Python at ordinary, reverse-flow and contact trial states.
 */
#include "kernels.h"
#include <math.h>
#include <stddef.h>

#define RU 8.31446261815324
#define MOLAR_MASS 0.004002602
#define TC 5.1953
#define PC 227460.0
#define OMEGA (-0.382)
static const double PI=3.1415926535897932384626433832795;
static const double pr_a=.457235583*RU*RU*TC*TC/PC;
static const double pr_b=.07779607*RU*TC/PC;
static const double kappa=.37464+1.54226*OMEGA-.26992*OMEGA*OMEGA;
static const double rg=RU/MOLAR_MASS;

static double cube_root(double x) { return x==0 ? 0 : copysign(pow(fabs(x),1.0/3.0),x); }
static void attraction(double T, double *a, double *da, double *dda) {
    double tr=T/TC, sr=sqrt(tr), base=1+kappa*(1-sr);
    *a=pr_a*base*base;
    *da=pr_a*(-base*kappa/(TC*sr));
    *dda=pr_a*kappa/(2*TC*TC)*(kappa/tr+base/(tr*sr));
}
static double z_factor(double p,double T) {
    double a,da,dda; attraction(T,&a,&da,&dda);
    double A=a*p/(RU*RU*T*T), B=pr_b*p/(RU*T);
    double ca=-(1-B),cb=A-3*B*B-2*B,cc=-(A*B-B*B-B*B*B);
    double pp=cb-ca*ca/3,qq=2*ca*ca*ca/27-ca*cb/3+cc;
    double disc=pow(qq/2,2)+pow(pp/3,3),off=-ca/3,roots[3]; int n;
    if(disc>1e-14) { roots[0]=cube_root(-qq/2+sqrt(disc))+cube_root(-qq/2-sqrt(disc))+off;n=1; }
    else if(fabs(disc)<=1e-14) { double u=cube_root(-qq/2);roots[0]=2*u+off;roots[1]=-u+off;n=2; }
    else { if(pp>=0) return NAN; double radius=2*sqrt(-pp/3);
        double arg=(3*qq/(2*pp))*sqrt(-3/pp),theta=acos(fmax(-1,fmin(1,arg)))/3;
        for(int i=0;i<3;i++) { roots[i]=radius*cos(theta-2*PI*i/3)+off; }
        n=3;
    }
    double z=-INFINITY;for(int i=0;i<n;i++) if(roots[i]>B && isfinite(roots[i])) z=fmax(z,roots[i]);
    return isfinite(z)?z:NAN;
}
static double density(double p,double T) { return MOLAR_MASS/(z_factor(p,T)*RU*T/p); }
static double log_volume(double rho) {
    double v=MOLAR_MASS/rho,sq=sqrt(2.0);
    return log((v+(1+sq)*pr_b)/(v+(1-sq)*pr_b))/(2*sq*pr_b*MOLAR_MASS);
}
static double u_departure(double T,double rho) {
    double a,da,dda;attraction(T,&a,&da,&dda);return (T*da-a)*log_volume(rho);
}
static double h_departure(double p,double T) {
    double a,da,dda;attraction(T,&a,&da,&dda);
    double z=z_factor(p,T),B=pr_b*p/(RU*T),sq=sqrt(2.0);
    double dep=RU*T*(z-1)+(T*da-a)*log((z+(1+sq)*B)/(z+(1-sq)*B))/(2*sq*pr_b);
    return dep/MOLAR_MASS;
}
static double u_ideal(double T) { return rg*(1.5*T-745.375); }
static double h_ideal(double T) { return rg*(2.5*T-745.375); }
static double temperature_u(double u) { return (u/rg+745.375)/1.5; }
static double temperature_h(double h) { return (h/rg+745.375)/2.5; }
static double pressure_rho(double T,double rho) {
    double v=MOLAR_MASS/rho,a,da,dda;attraction(T,&a,&da,&dda);
    if(v<=pr_b || T<=0) return NAN;
    return RU*T/(v-pr_b)-a/(v*(v+pr_b)+pr_b*(v-pr_b));
}
static NativeGas gas_properties(double m,double U,double V) {
    NativeGas g;g.rho=m/V;g.u=U/m;
    double T=fmax(temperature_u(g.u),2.2);
    for(int i=0;i<16;i++) { double next=fmax(temperature_u(g.u-u_departure(T,g.rho)),2.2);
        int done=fabs(next-T)<=1e-10*fmax(T,1);T=next;if(done) break; }
    g.T=T;g.p=pressure_rho(T,g.rho);g.h=h_ideal(T)+h_departure(g.p,T);return g;
}
static double temperature_ph(double p,double h) {
    double T=fmax(temperature_h(h),2.2);
    for(int i=0;i<16;i++) { double next=fmax(temperature_h(h-h_departure(p,T)),2.2);
        int done=fabs(next-T)<=1e-10*fmax(T,1);T=next;if(done) break; }return T;
}
static void local_isentropic(double p,double T,double *factor,double *exponent) {
    double rho=density(p,T),v=MOLAR_MASS/rho,a,da,dda;attraction(T,&a,&da,&dda);
    double d=v*(v+pr_b)+pr_b*(v-pr_b);
    double dpT=RU/(v-pr_b)-da/d;
    double dpR=(-RU*T/pow(v-pr_b,2)+a*2*(v+pr_b)/(d*d))*(-MOLAR_MASS/(rho*rho));
    double cv=1.5*rg+T*dda*log_volume(rho);
    double cp=cv+T*dpT*dpT/(rho*rho*dpR),gamma=cp/cv;
    *factor=p/(rho*dpR*gamma);*exponent=p*(gamma-1)/(gamma*T*dpT);
}
static double isentropic(double p,double T,double pd) {
    double f,e,fd,ed;local_isentropic(p,T,&f,&e);if(pd>=p) return f;
    double Td=fmax(T*pow(fmax(pd/p,1e-12),e),2.2);
    local_isentropic(fmax(pd,1),Td,&fd,&ed);return .5*(f+fd);
}
static double subsonic_cm(double r,double gamma,double rho,double T,double p) {
    return sqrt(fmax(2/(1-gamma)*rho*T/p*(pow(r,2*gamma)-pow(r,1+gamma)),0));
}
static void valve(double p,double pd,double T,double *cm,double *vel) {
    p=fmax(p,1);pd=fmax(fmin(pd,p),0);T=fmax(T,1);
    double g=fmax(1e-9,fmin(1-1e-9,isentropic(p,T,pd))),rho=fmax(density(p,T),1e-12);
    double r=fmax(pd/p,0),critical=pow(2*g/(g+1),1/(1-g)),eff;
    if(r<=critical) { eff=critical;*cm=sqrt(2/(1+g)*rho*T/p)*pow(2*g/(g+1),g/(1-g));*vel=sqrt(2/(1+g)*p/rho); }
    else { eff=r;*cm=subsonic_cm(r,g,rho,T,p);*vel=sqrt(fmax(2/(1-g)*p/rho*(1-pow(r,1-g)),0)); }
    double ref=subsonic_cm(.9999,g,rho,T,p);
    if(*cm>0 && ref>0) { double smooth=tanh(fmax(12*fabs(*cm/ref)*log(eff)/log(.9999),0));*cm*=smooth;*vel*=smooth; }
}

int native_gas_init(double p, double T, double volume, double *mU) {
    if (!(p > 0 && T >= 2.2 && volume > 0)) return 0;
    double rho = density(p, T);
    mU[0] = rho * volume;
    mU[1] = mU[0] * (u_ideal(T) + u_departure(T, rho));
    return isfinite(mU[0]) && isfinite(mU[1]) && mU[0] > 0;
}

int native_gas(double m, double U, double volume, NativeGas *gas) {
    if (!(m > 0 && volume > 0) || !isfinite(U)) return 0;
    *gas = gas_properties(m, U, volume);
    return gas->p > 0 && isfinite(gas->p) && isfinite(gas->h);
}

int native_orifice(double p1, double p2, double h1, double h2,
                   double cq_area, double opening, double *flow,
                   double *cm, double *velocity) {
    int forward = p1 >= p2;
    double p = forward ? p1 : p2, pd = forward ? p2 : p1;
    double T = temperature_ph(fmax(p, 1), forward ? h1 : h2);
    valve(p, pd, T, cm, velocity);
    double sign = forward ? 1 : -1;
    *velocity *= sign;
    *flow = opening == 0 || fabs(p1-p2) <= 1e-8 ? 0 :
        sign * cq_area * opening * fmax(p, 1) * *cm / sqrt(fmax(T, 1));
    if (opening == 0) *velocity = 0;
    return isfinite(*flow) && isfinite(*cm) && isfinite(*velocity);
}

double native_contact(double penetration, double velocity, double stiffness,
                      double damping, double pdis, int signed_force) {
    if (penetration <= 0) return 0;
    double fraction = pdis > 0 ? -expm1(-penetration / pdis) : 1;
    double force = stiffness * penetration + fraction * damping * velocity;
    return signed_force == 1 ? force : fmax(force, 0);
}

void native_stop_motion(double x, double v, double lower, double upper,
                        double *acceleration, double *velocity) {
    double vt = 1e-12 * fmax(fabs(v), 1);
    if ((x <= lower + 1e-12*fmax(fabs(lower),1) && v <= vt && *acceleration <= 0) ||
        (x >= upper - 1e-12*fmax(fabs(upper),1) && v >= -vt && *acceleration >= 0)) {
        *acceleration = 0;
        *velocity = 0;
    }
}

double native_signal(double t, double start, int stages, int cyclic, const double *data) {
    double elapsed = fmax(t-start,0), duration = 0, offset = 0;
    for (int i=0;i<stages;i++) duration += data[16+i];
    if (cyclic && duration > 0) elapsed = fmod(elapsed,duration);
    for (int i=0;i<stages;i++) {
        double d = data[16+i];
        if (elapsed < offset+d || i == stages-1)
            return d <= 0 ? data[8+i] : data[i] + (elapsed-offset)/d*(data[8+i]-data[i]);
        offset += d;
    }
    return data[8+stages-1];
}

double native_signal_break(double t, double end, double start, int stages,
                           int cyclic, const double *data) {
    double duration=0, offset=0, result=end;
    for (int i=0;i<stages;i++) duration += data[16+i];
    for (int i=0;i<stages;i++) {
        double event=start+offset;
        if (cyclic && duration > 0 && event <= t) {
            double cycle=fmax(0,floor((t-event)/duration)+1);
            event += cycle*duration;
            if (event <= t) event += duration;
        }
        if (event > t) result=fmin(result,event);
        offset += data[16+i];
    }
    return result;
}

static double ideal_temperature(const NativeMedium *m, double energy, double c) {
    double delta = energy - c*m->Tref;
    if (fabs(m->slope) <= 1e-15) return m->Tref + delta/c;
    double root = sqrt(fmax(c*c + 2*m->slope*delta, 0));
    double a = (-c+root)/m->slope, b = (-c-root)/m->slope;
    return m->Tref + (fabs(a)<=fabs(b) ? a : b);
}
int native_medium_init(const NativeMedium *medium, double p, double T, double V,
                       int legacy_ideal_initial, double *mU) {
    if (!(p>0 && T>0 && V>0)) return 0;
    if (medium->real_helium && !legacy_ideal_initial) return native_gas_init(p,T,V,mU);
    double mass=p*V/(medium->R*T),dt=T-medium->Tref;
    double u=medium->real_helium?u_ideal(T):(medium->cp-medium->R)*T+.5*medium->slope*dt*dt;
    mU[0]=mass;mU[1]=mass*u;
    return isfinite(mU[0]) && isfinite(mU[1]);
}
double native_density(const NativeMedium *m, double p, double T) {
    return m->real_helium ? density(p,T) : p/(m->R*T);
}
double native_temperature_ph(const NativeMedium *m, double p, double h) {
    return m->real_helium ? temperature_ph(p,h) : ideal_temperature(m,h,m->cp);
}
double native_viscosity(const NativeMedium *m, double T, int diagnostic) {
    /* Retain the ABI argument; flow and diagnostics use the same property. */
    (void)diagnostic;
    if (m->real_helium)
        return 1e-7*exp(.7501594*log(T)+35.76324/T-2212.129/(T*T)+.9212635);
    return m->mu*pow(T/m->muT,1.5)*(m->muT+m->S)/(T+m->S);
}
int native_medium_gas(const NativeMedium *medium, double m, double U, double V, NativeGas *g) {
    if (medium->real_helium) return native_gas(m,U,V,g);
    if (!(m>0 && V>0)) return 0;
    g->u=U/m; g->rho=m/V;
    g->T=ideal_temperature(medium,g->u,medium->cp-medium->R);
    g->p=g->rho*medium->R*g->T;
    double dt=g->T-medium->Tref;
    g->h=medium->cp*g->T+.5*medium->slope*dt*dt;
    return g->T>0 && isfinite(g->p) && isfinite(g->h);
}
static void medium_valve(const NativeMedium *m, double p, double pd, double T, double *cm, double *vel) {
    if (m->real_helium) { valve(p,pd,T,cm,vel); return; }
    p=fmax(p,1); pd=fmax(fmin(pd,p),0); T=fmax(T,1);
    double cp=m->cp+m->slope*(T-m->Tref);
    double g=fmax(1e-9,fmin(1-1e-9,(cp-m->R)/cp)),rho=fmax(native_density(m,p,T),1e-12);
    double r=pd/p,critical=pow(2*g/(g+1),1/(1-g)),eff;
    if(r<=critical) { eff=critical;*cm=sqrt(2/(1+g)*rho*T/p)*pow(2*g/(g+1),g/(1-g));*vel=sqrt(2/(1+g)*p/rho); }
    else { eff=r;*cm=subsonic_cm(r,g,rho,T,p);*vel=sqrt(fmax(2/(1-g)*p/rho*(1-pow(r,1-g)),0)); }
    double ref=subsonic_cm(.9999,g,rho,T,p);
    if(*cm>0 && ref>0) { double smooth=tanh(fmax(12*fabs(*cm/ref)*log(eff)/log(.9999),0));*cm*=smooth;*vel*=smooth; }
}
int native_medium_orifice(const NativeMedium *m, double p1, double p2, double h1, double h2,
                          double area, double opening, double *q, double *cm, double *v) {
    double p=fmax(p1,p2),pd=fmin(p1,p2),sign=p1>=p2?1:-1;
    double T=fmax(native_temperature_ph(m,fmax(p,1),p1>=p2?h1:h2),1);
    medium_valve(m,p,pd,T,cm,v);
    *q=fabs(p1-p2)<=1e-8?0:sign*area*opening*fmax(p,1)*(*cm)/sqrt(T);
    *v*=fabs(opening)<=1e-12?0:sign;
    return isfinite(*q) && isfinite(*cm) && isfinite(*v);
}
static double pipe_rough_limit(double rr) {
    return rr>0 ? 1/pow(-2*log10(rr/3.7),2) : 0;
}
static double pipe_friction_prepared(double re, double rr, double rough_limit) {
    if(re<=0) return 64000000;
    double lam=64/re;
    if(re<=89.96829989) return lam;
    double smooth=1/pow(-1.8*log10(6.9/re),2),turb=smooth;
    if(rr>0) { double r=re*rr,weight=r*r/(r*r+180*180);turb+=weight*(rough_limit-smooth); }
    double trans=pow((re-89.96829989)/2741.96700831,8.37293695);
    return lam+trans/(1+trans)*(turb-lam);
}
static double pipe_friction(double re, double rr) {
    return pipe_friction_prepared(re,rr,pipe_rough_limit(rr));
}
double native_pipe_flow(const NativeMedium *m, double p1, double p2, double T,
                        double d, double length, double rr, int kind) {
    if(fabs(p1-p2)<=1e-8) return 0;
    double p=fmax(fmax(p1,p2),1),pd=fmin(p1,p2),sign=p1>p2?1:-1;
    T=fmax(T,1);
    double area=PI*d*d/4,mu=native_viscosity(m,T,0),den=PI*d*mu;
    /* PNL0001/2/3 share compressible flow and its near-equilibrium smoothing.
       PNL0003 differs in storage placement, not in the resistance law. */
    double cm,vel;medium_valve(m,p,pd,T,&cm,&vel);
    if(kind==0) {
        double lam=pow(area*p*cm,2)/(16*PI*mu*length*T);
        if(4*lam/den<=1000) return sign*lam;
    }
    double rough_limit=pipe_rough_limit(rr);
    double base=area*p*cm/sqrt(T),q=sqrt(d/(length*.02))*base;
    for(int i=0;i<(kind==0?64:16);i++) {
        double next=sqrt(d/(length*pipe_friction_prepared(4*fabs(q)/den,rr,rough_limit)))*base;
        if(fabs(next-q)<=fmax(1e-12,fabs(q)*1e-9)) return sign*next;
        q=.5*(q+next);
    }
    return sign*q;
}
double native_pipe_flow_cached(NativePipeCache *cache, const NativeMedium *m,
                               double p1, double p2, double T, double d,
                               double length, double rr, int kind) {
    if(cache->valid && cache->p1==p1 && cache->p2==p2 && cache->T==T &&
       cache->diameter==d && cache->length==length && cache->roughness==rr && cache->kind==kind &&
       cache->medium.real_helium==m->real_helium && cache->medium.R==m->R &&
       cache->medium.cp==m->cp && cache->medium.Tref==m->Tref && cache->medium.slope==m->slope &&
       cache->medium.mu==m->mu && cache->medium.muT==m->muT && cache->medium.S==m->S)
        return cache->flow;
    double result=native_pipe_flow(m,p1,p2,T,d,length,rr,kind);
    cache->valid=0;
    if(isfinite(result)) {
        cache->medium=*m;cache->p1=p1;cache->p2=p2;cache->T=T;
        cache->diameter=d;cache->length=length;cache->roughness=rr;cache->kind=kind;
        cache->flow=result;cache->valid=1;
    }
    return result;
}
void native_pipe_diagnostics(const NativeMedium *m, double q, double p, double T,
                             double d, double length, double rr, int diagnostic, double *r) {
    double area=PI*d*d/4,re=4*fabs(q)/(PI*d*native_viscosity(m,T,0)),ff=pipe_friction(re,rr);
    r[0]=4*fabs(q)/(PI*d*native_viscosity(m,T,diagnostic));
    r[1]=fabs(q)*sqrt(T)/fmax(sqrt(d/(length*ff))*area*p,1e-18);
    r[2]=q/(fmax(native_density(m,p,T),1e-12)*area);r[3]=ff;
}
double native_limit_force(double penetration, double velocity, double stiffness,
                          double damping, double depth, int signed_force) {
    if(penetration<=0) return 0;
    double force=stiffness*penetration+(depth>0?fmin(penetration/depth,1):1)*damping*velocity;
    return signed_force==1?force:fmax(force,0);
}