#include "component_properties_internal.h"
#include "component_constants_internal.h"
#include <math.h>

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));
}
/* Derivative of the existing friction blend, local to this scalar solve. */
static double pipe_friction_derivative(double re,double rr,double rough,double *df) {
    double lam=64/re,dl=-lam/re;
    if(re<=89.96829989){*df=dl;return lam;}
    double a=-1.8*log10(6.9/re),smooth=1/(a*a);
    double ds=-2*smooth/a*1.8/(log(10.0)*re),turb=smooth,dt=ds;
    if(rr>0){double r=re*rr,w=r*r/(r*r+180*180),dw=2*w*(1-w)/re;
        turb+=w*(rough-smooth);dt=ds*(1-w)+dw*(rough-smooth);}
    double z=pow((re-89.96829989)/2741.96700831,8.37293695),b=z/(1+z);
    double db=8.37293695*b*(1-b)/(re-89.96829989);
    *df=dl*(1-b)+b*dt+db*(turb-lam);return lam+b*(turb-lam);
}
static double pipe_checked_solution(double re,double f,double K,double rr,double rough,
                                    double flow_per_re,NativePipeSolve *status) {
    double target=sqrt(K/f),change=fabs(target-re);
    /* Separate the absolute and relative tests so a large flow scale cannot
       make both sides overflow and accidentally satisfy the stopping test. */
    if(!isfinite(target) || !isfinite(target*flow_per_re) ||
       !(change<=re*1e-10 || change*flow_per_re<=1e-13))return NAN;
    double value=target*target*pipe_friction_prepared(target,rr,rough);
    status->relative_residual=fabs(value/K-1);
    if(isfinite(value) && status->relative_residual<=1e-9) {
        status->converged=1;return target;
    }
    return NAN;
}
double native_pipe_resistance(double K,double rr,double flow_per_re,NativePipeSolve *status) {
    NativePipeSolve local={0,0,0,INFINITY};if(!status)status=&local;
    *status=local;
    if(!(K>=0 && rr>=0 && flow_per_re>0) || !isfinite(K) || !isfinite(rr) || !isfinite(flow_per_re))return NAN;
    if(K==0){status->converged=1;status->relative_residual=0;return 0;}
    if(K/64<=89.96829989) {
        double re=K/64;
        status->relative_residual=fabs(64*re/K-1);
        if(re>0 && isfinite(re*flow_per_re) && status->relative_residual<=1e-9) {
            status->converged=1;return re;
        }
        return NAN;
    }
    double rough=pipe_rough_limit(rr),lo=0,hi=fmax(sqrt(K/.02),1),df;
    if(!isfinite(rough) || !isfinite(hi))return NAN;
    int bracketed=0;
    for(int i=0;i<128;i++) {
        double value=hi*hi*pipe_friction_prepared(hi,rr,rough);
        if(!isfinite(value))return NAN;
        if(value>=K){bracketed=1;break;}
        lo=hi;hi*=2;
        if(!isfinite(hi))return NAN;
    }
    /* F(0)=-K; each positive lower endpoint was checked while expanding.
       Never start iteration with an unchecked or non-finite upper residual. */
    if(!bracketed)return NAN;
    double re=fmax(lo,fmin(sqrt(K/.02),hi)),previous_residual=INFINITY;
    int previous_newton=0;
    for(int i=0;i<128;i++) {
        status->iterations=i+1;
        double f=pipe_friction_derivative(re,rr,rough,&df),value=re*re*f,F=value-K;
        if(!(f>0) || !isfinite(f) || !isfinite(value))return NAN;
        status->relative_residual=fabs(value/K-1);
        double target=pipe_checked_solution(re,f,K,rr,rough,flow_per_re,status);
        if(status->converged)return target;
        if(F>0)hi=re;else lo=re;
        double slope=2*re*f+re*re*df;
        double next=slope>0 && isfinite(slope) ? re-F/slope : NAN;
        /* An in-bracket Newton step is useful only if it reduces |F|.
           After a step that fails to halve it, bisect the retained bracket;
           successive Newton steps therefore cannot stagnate near an endpoint. */
        int bisect=(previous_newton && fabs(F)>.5*previous_residual) ||
                   !(slope>0) || !isfinite(slope) || !isfinite(next) || next<=lo || next>=hi;
        if(bisect) {
            next=lo+.5*(hi-lo);status->bisections++;
            if(next<=lo || next>=hi) {
                /* Adjacent floating-point endpoints: inspect both actual
                   candidates, then fail if neither meets the original tests. */
                double endpoints[]={lo,hi};
                for(int j=0;j<2;j++) {
                    double endpoint=endpoints[j];
                    target=pipe_checked_solution(endpoint,pipe_friction_prepared(endpoint,rr,rough),
                                                 K,rr,rough,flow_per_re,status);
                    if(status->converged)return target;
                }
                return NAN;
            }
        }
        previous_newton=!bisect;previous_residual=fabs(F);
        re=next;
    }
    return NAN;
}
double native_pipe_flow(const NativeMedium *m, double p1, double p2, double T,
                        double d, double length, double rr, int kind) {
    return native_pipe_flow_context(NULL,m,p1,p2,T,d,length,rr,kind);
}
double native_pipe_flow_context(NativePropertyCache *cache,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);
    NativePropertyState scratch,*up=property_pt(cache,m,p,T,&scratch);
    double area=PI*d*d/4,mu=property_viscosity(up),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;state_valve(cache,up,pd,&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 base=area*p*cm/sqrt(T),K=pow(4*base/den,2)*d/length;
    return sign*native_pipe_resistance(K,rr,den/4,NULL)*den/4;
}
double native_pipe_flow_cached(NativePipeCache *cache, const NativeMedium *m,
                               double p1, double p2, double T, double d,
                               double length, double rr, int kind) {
    return native_pipe_flow_cached_context(NULL,cache,m,p1,p2,T,d,length,rr,kind);
}
double native_pipe_flow_cached_context(NativePropertyCache *properties,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 &&
       same_medium(&cache->medium,m))
        return cache->flow;
    double result=properties?native_pipe_flow_context(properties,m,p1,p2,T,d,length,rr,kind):
                             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) {
    native_pipe_diagnostics_context(NULL,m,q,p,T,d,length,rr,diagnostic,r);
}
void native_pipe_diagnostics_context(NativePropertyCache *cache,const NativeMedium *m,double q,double p,double T,
                                     double d,double length,double rr,int diagnostic,double *r) {
    (void)diagnostic;
    NativePropertyState scratch,*state=property_pt(cache,m,p,T,&scratch);
    double area=PI*d*d/4,re=4*fabs(q)/(PI*d*property_viscosity(state)),ff=pipe_friction(re,rr);
    r[0]=re;
    r[1]=fabs(q)*sqrt(T)/fmax(sqrt(d/(length*ff))*area*p,1e-18);
    r[2]=q/(fmax(property_density(state),1e-12)*area);r[3]=ff;
}