#include "sopnamsp.h" #include "machdefs.h" #include #include #include #include #include "nbtri.h" #include "cspline.h" /*! \class SOPHYA::CSpline \ingroup NTools Spline 3 smoother (interpolator) for 1-D data points (y = f(x)) \sa CSpline2 */ /*! Createur pour spline 3 sur "x[0->n],y[0->n]" avec "yp1,ypn" derivees au premier et dernier points et "natural" indiquant les types de contraintes sur les derivees 2sd au premier et dernier point. "order" doit etre mis a "true" si le tableau de "x[]" n'est pas ordonne dans l'ordre des "x" croissants ("x[i]3 ); // allocation des buffers Y2 et tmp if( n>Nel || force ) { if( Y2 != NULL ) delete [] Y2; if( tmp != NULL ) delete [] tmp; Y2 = new double[n]; tmp = new double[n]; } // des-allocation eventuelle de X,Y if( XY_Created ) { if( !order || n>Nel || force ) { if( X != NULL ) delete [] X; X = NULL; if( Y != NULL ) delete [] Y; Y = NULL; XY_Created = false; } } // allocation eventuelle de X,Y if( order ) { if( !XY_Created || n>Nel || force ) { X = new double[n]; Y = new double[n]; XY_Created = true; } } Nel = n; corrupt_Y2 = true; if( x==NULL || y==NULL ) return; // Classement eventuel par ordre des x croissants if( order ) { if( tmp == NULL ) tmp = new double[n]; ind = (int_4 *) tmp; tri_double(x,ind,(int_4) Nel); for(int i=0;i0 ) if( X[i-1]>= X[i] ) { printf("CSpline::SetNewTab_Erreur: X[%d]>=X[%d] (%g>=%g)\n" ,i-1,i,X[i-1],X[i]); throw ParmError(PExcLongMessage("")); } } } else { X = x; Y = y; } } ////////////////////////////////////////////////////////////////////////////// //! destruction des divers tableaux en tenant compte des allocations/connections void CSpline::DelTab() { if( X != NULL && XY_Created ) delete [] X; X = NULL; if( Y != NULL && XY_Created ) delete [] Y; Y = NULL; if( Y2 != NULL ) delete [] Y2; Y2 = NULL; if( tmp != NULL ) delete [] tmp; tmp = NULL; } ////////////////////////////////////////////////////////////////////////////// /*! Pour changer les valeurs des derivees 1ere au 1er et dernier points Valeurs imposees des derivees 1ere au points "X[0]" et "X[Nel-1]". */ void CSpline::SetBound1er(double yp1,double ypn) { if( yp1 == YP1 && ypn == YPn ) return; YP1 = yp1; YPn = ypn; corrupt_Y2 = true; } //! Pour calculer les tableaux de coeff permettant le calcul des interpolations spline. void CSpline::ComputeCSpline() { // on ne fait rien si les tableaux ne sont pas connectes if( X == NULL || Y == NULL ) { printf("CSpline::ComputeCSpline()_Erreur: tableaux non connectes X=%p Y=%p\n" ,X,Y); return; } // On ne fait rien si rien n'a change! if( ! corrupt_Y2 ) return; // protection si tmp a ete desalloue pour gain de place (ex: CSpline2) if( tmp == NULL ) tmp = new double[Nel]; double p,qn,sig,un; if (Natural & Natural1) Y2[0] = tmp[0] = 0.0; else { Y2[0] = -0.5; tmp[0] = (3.0/(X[1]-X[0]))*((Y[1]-Y[0])/(X[1]-X[0])-YP1); } for (int i=1;i=0;k--) Y2[k] = Y2[k]*Y2[k+1] + tmp[k]; corrupt_Y2 = false; } //! Interpolation spline en \b x double CSpline::CSplineInt(double x) { int klo,khi,k; double h,b,a,y = 0.; if( corrupt_Y2 ) { cout<<"CSpline::CSplineInt: calcul des coef du spline corrupted"< 1) { k = (khi+klo) >> 1; if (X[k] > x) khi=k; else klo=k; } h=X[khi]-X[klo]; if (h == 0.0) { cout<<"CSpline::CSplineInt: pout khi="<3 && n2>3 ); int n = ( n1 < n2 ) ? n2 : n1; // allocation des buffers Y2 et tmp et des CSpline 1D if( n1>Nel1 || n2>Nel2 || force ) { if( Y2 != NULL ) delete [] Y2; if( tmp != NULL ) delete [] tmp; Y2 = new double[n1*n2]; tmp = new double[n]; // et les CSpline[n1] pour memoriser les interpolations sur x1(0->n1) if( S != NULL ) { for(int i=0;iFree_Tmp(); } Nel_S = n2; if( S != NULL ) { delete Sint; Sint = NULL;} Sint = new CSpline(n2,NULL,NULL,0.,0.,Natural); } // des-allocation eventuelle de X1,X2,Y if( XY_Created ) { if( !order || n1>Nel1 || n2>Nel2 || force ) { if( X1 != NULL ) delete [] X1; X1 = NULL; if( X2 != NULL ) delete [] X2; X2 = NULL; if( Y != NULL ) delete [] Y; Y = NULL; XY_Created = false; } } // allocation eventuelle de X1,X2,Y if( order ) { if( !XY_Created || n1>Nel1 || n2>Nel1 || force ) { X1 = new double[n1]; X2 = new double[n2]; Y = new double[n1*n2]; XY_Created = true; } } Nel1 = n1; Nel2 = n2; corrupt_Y2 = true; if( x1==NULL || x2==NULL || y==NULL ) return; // Classement eventuel par ordre des x1 et x2 croissants if( order ) { if( tmp == NULL ) tmp = new double[n]; ind = (int_4 *) tmp; double* Ytmp = new double[n1*n2]; // tri par valeur croissantes de x1 tri_double(x1,ind,(int_4) Nel1); for(int i=0;i0 ) if( X1[i-1] >= X1[i] ) { printf("CSpline::SetNewTab_Erreur: X1[%d]>=X1[%d] (%g>=%g)\n" ,i-1,i,X1[i-1],X1[i]); throw ParmError(PExcLongMessage("")); } for(int j=0;j0 ) if( X2[j-1] >= X2[j] ) { printf("CSpline::SetNewTab_Erreur: X2[%d]>=X2[%d] (%g>=%g)\n" ,j-1,j,X2[j-1],X2[j]); throw ParmError(PExcLongMessage("")); } for(int i=0;itmp = tmp; // connection de X1,Y au spline 1D sans ordre demande S[j]->SetNewTab(Nel1,X1,&Y[j*Nel1],false,false); // calcul des coeff splien pour l'interpolation future S[j]->ComputeCSpline(); } corrupt_Y2 = false; } ////////////////////////////////////////////////////////////////////////////// //! Calcule la valeur interpole (spline) pour le point \b (x1,x2) double CSpline2::CSplineInt(double x1,double x2) { // calcul de la valeur Y pour x=x1 et remplissage du tampon tmp for(int j=0;jCSplineInt(x1); // connection X2,tmp pour interpolation selon x=x2 Sint->SetNewTab(Nel2,X2,tmp,false,false); // calcul des coeff pour interpolation selon X2 Sint->ComputeCSpline(); // Interpolation finale return Sint->CSplineInt(x2); }