orsa_coord.cc File Reference

#include "orsa_coord.h"
#include <cmath>
#include <iostream>

Include dependency graph for orsa_coord.cc:

Go to the source code of this file.

Namespaces

namespace  orsa

Functions

void Interpolate (const vector< VectorWithParameter > vx_in, const double x, Vector &v_out, Vector &err_v_out)


Function Documentation

void Interpolate const vector< VectorWithParameter >  vx_in,
const double  x,
Vector &  v_out,
Vector &  err_v_out
 

Definition at line 101 of file orsa_coord.cc.

00101                                                                                                               {
00102     
00103     // cerr << "Interpolate() call..." << endl;
00104     
00105     unsigned int i,j,m;
00106     
00107     unsigned int i_closest;
00108     
00109     double      diff,denom;
00110     double      ho,hp;
00111     double      tmp_double;
00112 
00113     unsigned int n_points = vx_in.size();
00114     
00115     /* cerr << " -----> n_points: " << n_points << endl;
00116        for (j=0;j<n_points;j++)
00117        std::cerr << "vx_in[" << j << "].par = " <<  vx_in[j].par << "\n";
00118     */
00119     
00120     if (n_points < 2) {
00121       cerr << "too few points..." << endl;
00122       for (j=0;j<n_points;j++) {
00123         v_out = vx_in[0];
00124         err_v_out.Set(0,0,0);
00125       }
00126       return;
00127     }
00128     
00129     Vector * c = new Vector[n_points];
00130     Vector * d = new Vector[n_points];
00131     
00132     Vector w;
00133     
00134     diff = fabs(x-vx_in[0].par);
00135     i_closest = 0;
00136     
00137     for (j=0;j<n_points;j++) { 
00138       
00139       tmp_double = fabs(x-vx_in[j].par);
00140       if (tmp_double <  diff) {
00141         diff = tmp_double;
00142         i_closest = j;
00143       }
00144       
00145       c[j] = vx_in[j];
00146       d[j] = vx_in[j]; 
00147     }
00148     
00149     v_out     = vx_in[i_closest];
00150     err_v_out = vx_in[i_closest];
00151     
00152     // cerr << "local..  v: " << Length(v_out) << "  i_closest: " << i_closest << endl;
00153     
00154     i_closest--;
00155     
00156     for (m=1;m<=(n_points-2);m++) {
00157       for (i=0;i<(n_points-m);i++) {
00158         ho = vx_in[i].par   - x;
00159         hp = vx_in[i+m].par - x;
00160 
00161         denom = ho - hp;
00162 
00163         if (denom == 0.0) {
00164           cerr << "interpolate() --> Error: divide by zero\n";
00165           cerr << "i: " << i << "   m: " << m << endl;
00166           for (j=0;j<n_points;j++)
00167             cerr << "vx_in[" << j << "].par = " <<  vx_in[j].par << "\n";
00168           exit(0);
00169         }
00170         
00171         w = c[i+1] - d[i];
00172 
00173         d[i] = hp*w/denom;
00174         c[i] = ho*w/denom;
00175 
00176         // cerr << " ++++++++++  d[" << i << "]: " << Length(d[i]) << endl;
00177         // cerr << " ++++++++++  c[" << i << "]: " << Length(c[i]) << endl;
00178       }
00179       
00180       if ( (2*i_closest) < (n_points-m) ) {
00181         err_v_out = c[i_closest+1]; 
00182       } else {
00183         err_v_out = d[i_closest];
00184         i_closest--;
00185       }
00186       
00187       v_out += err_v_out;    
00188     }
00189     
00190     delete [] c;
00191     delete [] d;
00192   }


Generated on Tue Jan 11 15:26:42 2005 for liborsa by  doxygen 1.4.0