Mathc initiation/002n

Installer ce fichier dans votre répertoire de travail.

	x_juvw.h

/* --------------------------------- */
/* save as x_juvw.h                  */
/* --------------------------------- */
double Ju_uvw(
double (*P_f)(double x, double y, double z),

double (*P_X)(double u, double v, double w),
double (*P_Y)(double u, double v, double w),
double (*P_Z)(double u, double v, double w),

double (*Pu0)(double v, double w),
double (*Pu1)(double v, double w),
   int   nu,
   
double v,
   
double w
)
{
pt3d   p;

double u = 0;
   int i = 0;
double m = 0.;
double M = 0.;

 for(i = 0; i <= nu; i++)
 {
       if(i ==0 || i== nu){m = 1.;}
  else if(fmod(i,2) ==  0){m = 2.;}
  else                    {m = 4.;}
  
  u = ((*Pu0)(v,w)) + i*(((*Pu1)(v,w))-((*Pu0)(v,w)))/nu;

  p.x = u;
  p.y = v;   
  p.z = w;  
  
  M += m * (*P_f)( ((*P_X)(u,v,w)), 
                   ((*P_Y)(u,v,w)), 
                   ((*P_Z)(u,v,w)))
                       *
                   /* Jacobian */   
                      
                    fabs(
                    
              fxyz_x((*P_X),H,p) 
                       * 
   ( fxyz_y((*P_Y),H,p)*fxyz_z((*P_Z),H,p)
                       -
     fxyz_z((*P_Y),H,p)*fxyz_y((*P_Z),H,p) ) 
     
                       - 
                                       
              fxyz_y((*P_X),H,p) 
                       * 
   ( fxyz_x((*P_Y),H,p)*fxyz_z((*P_Z),H,p)
                       -
     fxyz_z((*P_Y),H,p)*fxyz_x((*P_Z),H,p) )
     
                       + 
                                       
              fxyz_z((*P_X),H,p) 
                       * 
   ( fxyz_x((*P_Y),H,p)*fxyz_y((*P_Z),H,p)
                       -
     fxyz_y((*P_Y),H,p)*fxyz_x((*P_Z),H,p) )
                        );  
 }

  return( ((((*Pu1)(v,w)) -((*Pu0)(v,w)))*M) / (3*nu) );
}
/* --------------------------------- */
double Jv_uvw(
double (*P_f)(double x, double y, double z),

double (*P_X)(double u, double v, double w),
double (*P_Y)(double u, double v, double w),
double (*P_Z)(double u, double v, double w),
   
double (*Pu0)(double v, double w),
double (*Pu1)(double v, double w),
   int   nu,

double (*Pv0)(double u),
double (*Pv1)(double u),
   int   nv,
   
double w
)
{
   int i = 0;
double m = 0.;
double M = 0.;

 for(i = 0; i <= nv; i++)
 {
       if(i ==0 || i== nv){m = 1.;}
  else if(fmod(i,2) ==  0){m = 2.;}
  else                    {m = 4.;}

  M += m * Ju_uvw((*P_f),
                       (*P_X),(*P_Y),(*P_Z),   
                       (*Pu0),(*Pu1),nu,                               
                     (((*Pv0)(w))+i*(((*Pv1)(w))-((*Pv0)(w)))/nv),w);
 }

 return( ((((*Pv1)(w)) -((*Pv0)(w)))*M) / (3*nv) );
}
/* --------------------------------- */
double J_uvw(
double (*P_f)(double x, double y, double z),

double (*P_X)(double u, double v, double w),
double (*P_Y)(double u, double v, double w),
double (*P_Z)(double u, double v, double w),
   
double (*Pu0)(double v, double w),
double (*Pu1)(double v, double w),
   int    nu,
   
double (*Pv0)(double u),
double (*Pv1)(double u),
   int   nv,

double w0,
double w1,
   int nw
)
{
   int i = 0;
double m = 0.;
double M = 0.;

 for(i = 0; i <= nw; i++)
 {
       if(i ==0 || i== nw){m = 1.;}
  else if(fmod(i,2) ==  0){m = 2.;}
  else                    {m = 4.;}

  M += m *       Jv_uvw((*P_f),
                             (*P_X),(*P_Y),(*P_Z),  
                             (*Pu0),(*Pu1),nu,                             
                             (*Pv0),(*Pv1),nv,                                
                             (w0+i*(w1-w0)/nw));
 }

  return( ((w1 -w0)*M) / (3*nw) );
}
/* --------------------------------- */
/* --------------------------------- */