优化端口传递物性参数方式,减少上下游元件重复计算量

This commit is contained in:
ljz committed 2026-09-11 14:42:19 +08:00
1 parent 91bd9fb252
commit 5d5a2e1843
12 files changed
+773 -65

No files matched your search

+168 -42
View File
@@ -16,6 +16,71 @@ 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 const NativeMedium helium_medium={1,RU/MOLAR_MASS,2.5*RU/MOLAR_MASS,293.15,0,1.96e-5,293.15,79.4};
static int same_medium(const NativeMedium *a,const NativeMedium *b) {
return a->real_helium==b->real_helium && a->R==b->R && a->cp==b->cp &&
a->Tref==b->Tref && a->slope==b->slope && a->mu==b->mu && a->muT==b->muT && a->S==b->S;
}
void native_properties_init(NativePropertyCache *cache,NativePropertyState *states,size_t capacity) {
cache->states=states;cache->count=0;cache->capacity=states?capacity:0;
}
static NativePropertyState *property_new(NativePropertyCache *cache,const NativeMedium *m,
double p,double T,NativePropertyState *scratch) {
int valid=p>0 && T>0 && isfinite(p) && isfinite(T);
NativePropertyState *s=valid && cache && cache->count<cache->capacity ? &cache->states[cache->count++] : scratch;
*s=(NativePropertyState){0};s->medium=*m;s->p=p;s->T=T;s->valid=valid?NATIVE_PROPERTY_PT:0;
return s;
}
static NativePropertyState *property_pt(NativePropertyCache *cache,const NativeMedium *m,
double p,double T,NativePropertyState *scratch) {
if(cache)for(size_t i=0;i<cache->count;i++) {
NativePropertyState *s=&cache->states[i];
if((s->valid&NATIVE_PROPERTY_PT) && s->p==p && s->T==T && same_medium(&s->medium,m))return s;
}
return property_new(cache,m,p,T,scratch);
}
static double property_density(NativePropertyState *s) {
if(!(s->valid&NATIVE_PROPERTY_RHO)) {
s->rho=native_density(&s->medium,s->p,s->T);
if(s->rho>0 && isfinite(s->rho))s->valid|=NATIVE_PROPERTY_RHO;
}
return s->rho;
}
static double property_viscosity(NativePropertyState *s) {
if(!(s->valid&NATIVE_PROPERTY_MU)) {
s->mu=native_viscosity(&s->medium,s->T,0);
if(s->mu>0 && isfinite(s->mu))s->valid|=NATIVE_PROPERTY_MU;
}
return s->mu;
}
static void remember_gas(NativePropertyCache *cache,const NativeMedium *m,const NativeGas *gas) {
/* Below helium's critical temperature retain the existing vapor-root/PH
selection. A future phase-aware medium contract can carry that state. */
if(!cache || (m->real_helium && gas->T<=TC))return;
NativePropertyState scratch,*s=property_pt(cache,m,gas->p,gas->T,&scratch);
if(!(s->valid&NATIVE_PROPERTY_PT) || !isfinite(gas->h))return;
if((s->valid&NATIVE_PROPERTY_H) && s->h!=gas->h)
s=property_new(cache,m,gas->p,gas->T,&scratch);
s->h=gas->h;s->valid|=NATIVE_PROPERTY_H;
if(gas->rho>0 && isfinite(gas->rho)){s->rho=gas->rho;s->valid|=NATIVE_PROPERTY_RHO;}
}
double native_temperature_ph_context(NativePropertyCache *cache,const NativeMedium *m,double p,double h) {
if(cache)for(size_t i=0;i<cache->count;i++) {
NativePropertyState *s=&cache->states[i];
if((s->valid&NATIVE_PROPERTY_H) && s->p==p && s->h==h && same_medium(&s->medium,m))return s->T;
}
double T=native_temperature_ph(m,p,h);
if(cache && isfinite(h) && T>0 && isfinite(T) && p>0 && isfinite(p)) {
NativePropertyState scratch,*s=property_pt(cache,m,p,T,&scratch);
if((s->valid&NATIVE_PROPERTY_H) && s->h!=h)s=property_new(cache,m,p,T,&scratch);
s->h=h;s->valid|=NATIVE_PROPERTY_H;
}
return T;
}
double native_density_context(NativePropertyCache *cache,const NativeMedium *m,double p,double T) {
NativePropertyState scratch;return property_density(property_pt(cache,m,p,T,&scratch));
}
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) {
@@ -68,39 +133,52 @@ static NativeGas gas_properties(double m,double U,double V) {
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;
g.T=T;g.p=pressure_rho(T,g.rho);
/* Same energy reference as u_ideal/h_ideal. Avoid a redundant cubic solve
in the single-root region; preserve vapor-root semantics below TC. */
g.h=T>TC ? g.u+g.p/g.rho : 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);
static void local_isentropic(NativePropertyState *s) {
if(s->valid&NATIVE_PROPERTY_ISENTROPIC)return;
double p=s->p,T=s->T,rho=property_density(s),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);
s->isentropic_factor=p/(rho*dpR*gamma);s->isentropic_exponent=p*(gamma-1)/(gamma*T*dpT);
if(isfinite(s->isentropic_factor) && isfinite(s->isentropic_exponent))s->valid|=NATIVE_PROPERTY_ISENTROPIC;
}
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 isentropic(NativePropertyCache *cache,NativePropertyState *up,double pd) {
local_isentropic(up);if(pd>=up->p)return up->isentropic_factor;
double Td=fmax(up->T*pow(fmax(pd/up->p,1e-12),up->isentropic_exponent),2.2);
NativePropertyState scratch,*down=property_pt(cache,&up->medium,fmax(pd,1),Td,&scratch);
local_isentropic(down);return .5*(up->isentropic_factor+down->isentropic_factor);
}
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);
static void state_valve(NativePropertyCache *cache,NativePropertyState *up,double pd,double *cm,double *vel) {
double p=up->p,T=up->T;pd=fmax(fmin(pd,p),0);
const NativeMedium *m=&up->medium;
double cp=m->cp+m->slope*(T-m->Tref);
double factor=m->real_helium?isentropic(cache,up,pd):(cp-m->R)/cp;
double g=fmax(1e-9,fmin(1-1e-9,factor)),rho=fmax(property_density(up),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; }
}
static void medium_valve(NativePropertyCache *cache,const NativeMedium *m,double p,double pd,double T,double *cm,double *vel) {
NativePropertyState scratch,*up=property_pt(cache,m,fmax(p,1),fmax(T,1),&scratch);
state_valve(cache,up,pd,cm,vel);
}
int native_gas_init(double p, double T, double volume, double *mU) {
if (!(p > 0 && T >= 2.2 && volume > 0)) return 0;
@@ -115,14 +193,21 @@ int native_gas(double m, double U, double volume, NativeGas *gas) {
*gas = gas_properties(m, U, volume);
return gas->p > 0 && isfinite(gas->p) && isfinite(gas->h);
}
int native_gas_context(NativePropertyCache *cache,double m,double U,double V,NativeGas *gas) {
int ok=native_gas(m,U,V,gas);if(ok)remember_gas(cache,&helium_medium,gas);return ok;
}
int native_orifice(double p1, double p2, double h1, double h2,
double cq_area, double opening, double *flow,
double *cm, double *velocity) {
return native_orifice_context(NULL,p1,p2,h1,h2,cq_area,opening,flow,cm,velocity);
}
int native_orifice_context(NativePropertyCache *cache,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 T = native_temperature_ph_context(cache,&helium_medium,fmax(p,1),forward?h1:h2);
medium_valve(cache,&helium_medium,p,pd,T,cm,velocity);
double sign = forward ? 1 : -1;
*velocity *= sign;
*flow = opening == 0 || fabs(p1-p2) <= 1e-8 ? 0 :
@@ -218,22 +303,18 @@ int native_medium_gas(const NativeMedium *medium, double m, double U, double V,
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_gas_context(NativePropertyCache *cache,const NativeMedium *m,double mass,double U,double V,NativeGas *g) {
int ok=native_medium_gas(m,mass,U,V,g);if(ok)remember_gas(cache,m,g);return ok;
}
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) {
return native_medium_orifice_context(NULL,m,p1,p2,h1,h2,area,opening,q,cm,v);
}
int native_medium_orifice_context(NativePropertyCache *cache,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);
double T=fmax(native_temperature_ph_context(cache,m,fmax(p,1),p1>=p2?h1:h2),1);
medium_valve(cache,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);
@@ -253,38 +334,77 @@ static double pipe_friction_prepared(double re, double rr, double rough_limit) {
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);
}
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/64<=89.96829989){status->converged=1;status->relative_residual=0;return K/64;}
double rough=pipe_rough_limit(rr),lo=0,hi=fmax(sqrt(K/.02),1),df;
for(int i=0;i<128 && hi*hi*pipe_friction_prepared(hi,rr,rough)<K;i++)hi*=2;
if(!isfinite(hi))return NAN;
double re=fmin(sqrt(K/.02),hi);
for(int i=0;i<80;i++) {
status->iterations=i+1;
double f=pipe_friction_derivative(re,rr,rough,&df),F=re*re*f-K;
double target=sqrt(K/f);
if(fabs(target-re)*flow_per_re<=fmax(1e-13,re*flow_per_re*1e-10)) {
/* Check the returned point, not just the change of successive guesses. */
status->relative_residual=fabs(target*target*pipe_friction_prepared(target,rr,rough)/K-1);
if(isfinite(target) && status->relative_residual<=1e-9){status->converged=1;return target;}
}
if(F>0)hi=re;else lo=re;
double next=re-F/(2*re*f+re*re*df);
if(!isfinite(next) || next<=lo || next>=hi){next=.5*(lo+hi);status->bisections++;}
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);
double area=PI*d*d/4,mu=native_viscosity(m,T,0),den=PI*d*mu;
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;medium_valve(m,p,pd,T,&cm,&vel);
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 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 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 &&
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)
same_medium(&cache->medium,m))
return cache->flow;
double result=native_pipe_flow(m,p1,p2,T,d,length,rr,kind);
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;
@@ -295,10 +415,16 @@ double native_pipe_flow_cached(NativePipeCache *cache, const NativeMedium *m,
}
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));
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(native_density(m,p,T),1e-12)*area);r[3]=ff;
r[2]=q/(fmax(property_density(state),1e-12)*area);r[3]=ff;
}
double native_limit_force(double penetration, double velocity, double stiffness,
double damping, double depth, int signed_force) {
+35
View File
@@ -1,5 +1,6 @@
#ifndef NATIVE_KERNELS_H
#define NATIVE_KERNELS_H
#include <stddef.h>
typedef struct { double p, T, rho, u, h; } NativeGas;
typedef struct {
int velocity_index;
@@ -8,6 +9,30 @@ typedef struct {
} NativeStop;
/* Constants are emitted per medium instance by the model compiler. */
typedef struct { int real_helium; double R, cp, Tref, slope, mu, muT, S; } NativeMedium;
/* Caller-owned scratch for ONE model_eval. Keys use exact values and a copy of
* every medium constant. No process/thread global cache or approximate reuse. */
enum {
NATIVE_PROPERTY_PT=1, NATIVE_PROPERTY_H=2, NATIVE_PROPERTY_RHO=4,
NATIVE_PROPERTY_MU=8, NATIVE_PROPERTY_ISENTROPIC=16
};
typedef struct {
NativeMedium medium;
double p, T, h, rho, mu, isentropic_factor, isentropic_exponent;
unsigned valid;
} NativePropertyState;
typedef struct {
NativePropertyState *states;
size_t count, capacity;
} NativePropertyCache;
void native_properties_init(NativePropertyCache *, NativePropertyState *, size_t capacity);
int native_gas_context(NativePropertyCache *, double m, double U, double V, NativeGas *);
int native_medium_gas_context(NativePropertyCache *, const NativeMedium *, double m, double U, double V, NativeGas *);
double native_temperature_ph_context(NativePropertyCache *, const NativeMedium *, double p, double h);
double native_density_context(NativePropertyCache *, const NativeMedium *, double p, double T);
int native_orifice_context(NativePropertyCache *, double p1, double p2, double h1, double h2,
double area, double opening, double *q, double *cm, double *v);
int native_medium_orifice_context(NativePropertyCache *, const NativeMedium *, double p1, double p2, double h1, double h2,
double area, double opening, double *q, double *cm, double *v);
/* One entry per pipe branch, zero-initialized for each model_eval. Never shared
* across solver trials. Exact inputs, including medium constants, form the key. */
typedef struct {
@@ -15,6 +40,16 @@ typedef struct {
double p1, p2, T, diameter, length, roughness, flow;
int kind, valid;
} NativePipeCache;
typedef struct { int converged, iterations, bisections; double relative_residual; } NativePipeSolve;
/* Solve Re^2*f(Re)=K for the shared pipe resistance law. On failure returns
* NAN, with converged=0; flow_per_re converts the local stopping test to kg/s. */
double native_pipe_resistance(double K, double roughness, double flow_per_re, NativePipeSolve *status);
double native_pipe_flow_context(NativePropertyCache *, const NativeMedium *, double p1, double p2, double T,
double diameter, double length, double roughness, int kind);
double native_pipe_flow_cached_context(NativePropertyCache *, NativePipeCache *, const NativeMedium *,
double p1, double p2, double T, double diameter, double length, double roughness, int kind);
void native_pipe_diagnostics_context(NativePropertyCache *, const NativeMedium *, double q, double p, double T,
double diameter, double length, double roughness, int diagnostic, double *result);
int native_medium_init(const NativeMedium *, double p, double T, double V, int legacy_ideal_initial, double *mU);
double native_density(const NativeMedium *, double p, double T);
double native_temperature_ph(const NativeMedium *, double p, double h);