LORENE: mtbl_cf_val_point.C Source File

00001 /*
00002  * Member functions of the Mtbl_cf class for computing the value of a field
00003  *  at an arbitrary point
00004  *
00005  * (see file mtbl_cf.h for the documentation).
00006  */
00007 
00008 /*
00009  *   Copyright (c) 1999-2002 Eric Gourgoulhon
00010  *
00011  *   This file is part of LORENE.
00012  *
00013  *   LORENE is free software; you can redistribute it and/or modify
00014  *   it under the terms of the GNU General Public License as published by
00015  *   the Free Software Foundation; either version 2 of the License, or
00016  *   (at your option) any later version.
00017  *
00018  *   LORENE is distributed in the hope that it will be useful,
00019  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
00020  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00021  *   GNU General Public License for more details.
00022  *
00023  *   You should have received a copy of the GNU General Public License
00024  *   along with LORENE; if not, write to the Free Software
00025  *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00026  *
00027  */
00028 
00029 
00030 char mtbl_cf_val_point_C[] = "$Header: /cvsroot/Lorene/C++/Source/Mtbl_cf/mtbl_cf_val_point.C,v 1.13 2012/01/17 15:09:05 j_penner Exp $" ;
00031 
00032 /*
00033  * $Id: mtbl_cf_val_point.C,v 1.13 2012/01/17 15:09:05 j_penner Exp $
00034  * $Log: mtbl_cf_val_point.C,v $
00035  * Revision 1.13  2012/01/17 15:09:05  j_penner
00036  * using MAX_BASE_2 for the phi coordinate
00037  *
00038  * Revision 1.12  2009/10/08 16:21:16  j_novak
00039  * Addition of new bases T_COS and T_SIN.
00040  *
00041  * Revision 1.11  2007/12/20 09:11:08  jl_cornou
00042  * Correction of an error in op_sxpun about Jacobi(0,2) polynomials
00043  *
00044  * Revision 1.10  2007/12/11 15:28:16  jl_cornou
00045  * Jacobi(0,2) polynomials partially implemented
00046  *
00047  * Revision 1.9  2007/10/23 17:15:13  j_novak
00048  * Added the bases T_COSSIN_C and T_COSSIN_S
00049  *
00050  * Revision 1.8  2006/06/06 14:57:01  j_novak
00051  * Summation functions for angular coefficients at xi=+/-1.
00052  *
00053  * Revision 1.7  2006/05/30 08:30:15  n_vasset
00054  * Implementation of sine-like bases (T_SIN_P, T_SIN_I, T_COSSIN_SI, etc...).
00055  *
00056  * Revision 1.6  2005/05/27 14:55:00  j_novak
00057  * Added new bases T_COSSIN_CI and T_COS_I
00058  *
00059  * Revision 1.5  2005/02/16 15:10:39  m_forot
00060  * Correct the case T_COSSIN_C
00061  *
00062  * Revision 1.4  2004/12/17 13:35:03  m_forot
00063  * Add the case T_LEG
00064  *
00065  * Revision 1.3  2002/05/11 12:37:31  e_gourgoulhon
00066  * Added basis T_COSSIN_SI.
00067  *
00068  * Revision 1.2  2002/05/05 16:22:33  e_gourgoulhon
00069  * Added the case of the theta basis T_COSSIN_SP.
00070  *
00071  * Revision 1.1.1.1  2001/11/20 15:19:27  e_gourgoulhon
00072  * LORENE
00073  *
00074  * Revision 1.7  2000/09/08  16:26:43  eric
00075  * Ajout de la base T_SIN_I.
00076  *
00077  * Revision 1.6  2000/09/08  16:07:16  eric
00078  * Ajout de la base P_COSSIN_I
00079  *
00080  * Revision 1.5  2000/09/06  14:00:19  eric
00081  * Ajout de la base T_COS_I.
00082  *
00083  * Revision 1.4  1999/12/29  13:11:42  eric
00084  * Ajout de la fonction val_point_jk.
00085  *
00086  * Revision 1.3  1999/12/07  15:10:45  eric
00087  * *** empty log message ***
00088  *
00089  * Revision 1.2  1999/12/07  14:52:34  eric
00090  * Changement ordre des arguments (phi,theta,xi) --> (xi,theta,phi)
00091  *
00092  * Revision 1.1  1999/12/06  16:47:39  eric
00093  * Initial revision
00094  *
00095  *
00096  * $Header: /cvsroot/Lorene/C++/Source/Mtbl_cf/mtbl_cf_val_point.C,v 1.13 2012/01/17 15:09:05 j_penner Exp $
00097  *
00098  */
00099 
00100 // Headers Lorene
00101 #include "mtbl_cf.h"
00102 #include "proto.h"
00103 
00104     //-------------------------------------------------------------//
00105     //      version for an arbitrary point in (xi,theta',phi')     //
00106     //-------------------------------------------------------------//
00107 
00108 double Mtbl_cf::val_point(int l, double x, double theta, double phi) const {
00109 
00110 // Routines de sommation
00111 static void (*som_r[MAX_BASE])
00112         (double*, const int, const int, const int, const double, double*) ;
00113 static void (*som_tet[MAX_BASE])
00114         (double*, const int, const int, const double, double*) ;
00115 static void (*som_phi[MAX_BASE_2])
00116         (double*, const int, const double, double*) ;
00117 static int premier_appel = 1 ;
00118 
00119 // Initialisations au premier appel
00120 // --------------------------------
00121     if (premier_appel == 1) {
00122 
00123     premier_appel = 0 ;
00124 
00125     for (int i=0 ; i<MAX_BASE ; i++) {
00126         if(i%2 == 0){
00127         som_phi[i/2] = som_phi_pas_prevu ;
00128         }
00129         som_tet[i] = som_tet_pas_prevu ;
00130         som_r[i]   = som_r_pas_prevu ;
00131     }
00132 
00133     som_r[R_CHEB >> TRA_R] = som_r_cheb ;
00134     som_r[R_CHEBP >> TRA_R] = som_r_chebp ;
00135     som_r[R_CHEBI >> TRA_R] = som_r_chebi ;
00136     som_r[R_CHEBU >> TRA_R] = som_r_chebu ;
00137     som_r[R_CHEBPIM_P >> TRA_R] = som_r_chebpim_p ;
00138     som_r[R_CHEBPIM_I >> TRA_R] = som_r_chebpim_i ;
00139     som_r[R_CHEBPI_P >> TRA_R] = som_r_chebpi_p ;
00140     som_r[R_CHEBPI_I >> TRA_R] = som_r_chebpi_i ;
00141     som_r[R_JACO02 >> TRA_R] = som_r_jaco02 ;
00142 
00143     som_tet[T_COS >> TRA_T] = som_tet_cos ;
00144     som_tet[T_SIN >> TRA_T] = som_tet_sin ;
00145     som_tet[T_COS_P >> TRA_T] = som_tet_cos_p ;
00146     som_tet[T_COS_I >> TRA_T] = som_tet_cos_i ;
00147     som_tet[T_SIN_P >> TRA_T] = som_tet_sin_p ;
00148     som_tet[T_SIN_I >> TRA_T] = som_tet_sin_i ;
00149     som_tet[T_COSSIN_CP >> TRA_T] = som_tet_cossin_cp ;
00150     som_tet[T_COSSIN_CI >> TRA_T] = som_tet_cossin_ci ;
00151     som_tet[T_COSSIN_SP >> TRA_T] = som_tet_cossin_sp ;
00152     som_tet[T_COSSIN_SI >> TRA_T] = som_tet_cossin_si ;
00153     som_tet[T_COSSIN_C >> TRA_T] = som_tet_cossin_c ;
00154     som_tet[T_COSSIN_S >> TRA_T] = som_tet_cossin_s ;
00155 
00156     som_phi[P_COSSIN >> TRA_P] = som_phi_cossin ;
00157     som_phi[P_COSSIN_P >> TRA_P] = som_phi_cossin_p ;
00158     som_phi[P_COSSIN_I >> TRA_P] = som_phi_cossin_i ;
00159 
00160     }   // fin des operations de premier appel
00161 
00162 
00163     assert (etat != ETATNONDEF) ; 
00164 
00165     double resu ;           // valeur de retour    
00166 
00167 // Cas ou tous les coefficients sont nuls :
00168     if (etat == ETATZERO ) {
00169     resu = 0 ;
00170     return resu ; 
00171     } 
00172  
00173 // Nombre de points en phi, theta et r :
00174     int np = mg->get_np(l) ;     
00175     int nt = mg->get_nt(l) ;
00176     int nr = mg->get_nr(l) ;
00177 
00178 // Bases de developpement : 
00179     int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
00180     int base_t = (base.b[l] & MSQ_T) >> TRA_T ;
00181     int base_p = (base.b[l] & MSQ_P) >> TRA_P ;
00182     
00183 //--------------------------------------
00184 // Calcul de la valeur au point demande
00185 //--------------------------------------
00186 
00187 // Pointeur sur le tableau contenant les coefficients: 
00188 
00189     assert(etat == ETATQCQ) ; 
00190     Tbl* tbcf = t[l] ; 
00191     
00192     if (tbcf->get_etat() == ETATZERO ) {
00193     resu = 0 ;
00194     return resu ; 
00195     } 
00196 
00197 
00198     assert(tbcf->get_etat() == ETATQCQ) ; 
00199 
00200     double* cf = tbcf->t ;
00201 
00202     // Tableaux de travail 
00203     double* trp = new double [np+2] ;
00204     double* trtp = new double [(np+2)*nt] ;
00205  
00206     if (nr == 1) {
00207     
00208 // Cas particulier nr = 1 (Fonction purement angulaire) : 
00209 // ----------------------
00210     
00211     som_tet[base_t](cf, nt, np, theta, trp) ;   // sommation sur theta
00212     }
00213     else {
00214 
00215 // Cas general
00216 // -----------
00217 
00218     som_r[base_r](cf, nr, nt, np, x, trtp) ;    // sommation sur r
00219     som_tet[base_t](trtp, nt, np, theta, trp) ; // sommation sur theta
00220     }
00221 
00222 // Sommation sur phi
00223 // -----------------
00224 
00225     if (np == 1) {
00226     resu = trp[0] ;     // cas axisymetrique
00227     }
00228     else {
00229     som_phi[base_p](trp, np, phi, &resu) ;      // sommation sur phi
00230     }
00231 
00232     // Menage
00233     delete [] trp ;
00234     delete [] trtp ;
00235     
00236     // Termine
00237     return resu ;  
00238  
00239 }
00240  
00241 
00242 
00243     //-------------------------------------------------------------//
00244     //      version for an arbitrary point in xi           //
00245     //      but collocation point in (theta',phi')         //
00246     //-------------------------------------------------------------//
00247 
00248 double Mtbl_cf::val_point_jk(int l, double x, int j0, int k0) const {
00249 
00250 // Routines de sommation
00251 static void (*som_r[MAX_BASE])
00252         (double*, const int, const int, const int, const double, double*) ;
00253 static int premier_appel = 1 ;
00254 
00255 // Initialisations au premier appel
00256 // --------------------------------
00257     if (premier_appel == 1) {
00258 
00259     premier_appel = 0 ;
00260 
00261     for (int i=0 ; i<MAX_BASE ; i++) {
00262         som_r[i]   = som_r_pas_prevu ;
00263     }
00264 
00265     som_r[R_CHEB >> TRA_R] = som_r_cheb ;
00266     som_r[R_CHEBP >> TRA_R] = som_r_chebp ;
00267     som_r[R_CHEBI >> TRA_R] = som_r_chebi ;
00268     som_r[R_CHEBU >> TRA_R] = som_r_chebu ;
00269     som_r[R_CHEBPIM_P >> TRA_R] = som_r_chebpim_p ;
00270     som_r[R_CHEBPIM_I >> TRA_R] = som_r_chebpim_i ;
00271     som_r[R_CHEBPI_P >> TRA_R] = som_r_chebpi_p ;
00272     som_r[R_CHEBPI_I >> TRA_R] = som_r_chebpi_i ;
00273     som_r[R_JACO02 >> TRA_R] = som_r_jaco02 ;
00274 
00275     }   // fin des operations de premier appel
00276 
00277     assert (etat != ETATNONDEF) ; 
00278 
00279     double resu ;           // valeur de retour    
00280 
00281 // Cas ou tous les coefficients sont nuls :
00282     if (etat == ETATZERO ) {
00283     resu = 0 ;
00284     return resu ; 
00285     } 
00286  
00287 // Nombre de points en phi, theta et r :
00288     int np = mg->get_np(l) ;     
00289     int nt = mg->get_nt(l) ;
00290     int nr = mg->get_nr(l) ;
00291 
00292 // Bases de developpement : 
00293     int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
00294 
00295 //------------------------------------------------------------------------
00296 //  Valeurs des fonctions de base en phi aux points de collocation en phi
00297 //   et des fonctions de base en theta aux points de collocation en theta
00298 //------------------------------------------------------------------------
00299 
00300     Tbl tab_phi = base.phi_functions(l, np) ; 
00301     Tbl tab_theta = base.theta_functions(l, nt) ; 
00302 
00303     
00304 //--------------------------------------
00305 // Calcul de la valeur au point demande
00306 //--------------------------------------
00307 
00308 // Pointeur sur le tableau contenant les coefficients: 
00309 
00310     assert(etat == ETATQCQ) ; 
00311     Tbl* tbcf = t[l] ; 
00312     
00313     if (tbcf->get_etat() == ETATZERO ) {
00314     resu = 0 ;
00315     return resu ; 
00316     } 
00317 
00318 
00319     assert(tbcf->get_etat() == ETATQCQ) ; 
00320 
00321     double* cf = tbcf->t ;
00322 
00323     // Tableau de travail 
00324     double* coef_tp = new double [(np+2)*nt] ;
00325  
00326 
00327 // 1/ Sommation sur r
00328 // ------------------
00329 
00330     som_r[base_r](cf, nr, nt, np, x, coef_tp) ;
00331 
00332 
00333 // 2/ Sommation sur theta et phi
00334 // -----------------------------
00335     double* pi = coef_tp ;  // pointeur courant sur les coef en theta et phi
00336 
00337 // Sommation sur le premier phi, k=0
00338 
00339     double somt = 0 ;
00340     for (int j=0 ; j<nt ; j++) {
00341     somt += (*pi) * tab_theta(0, j, j0) ;
00342     pi++ ;  // theta suivant
00343     }
00344     resu = somt * tab_phi(0, k0)    ;
00345  
00346     if (np > 1) {   // sommation sur phi
00347 
00348     // On saute le phi suivant (sin(0)), k=1
00349     pi += nt ;
00350 
00351     // Sommation sur le reste des phi (pour k=2,...,np)
00352 
00353     int base_t = base.b[l] & MSQ_T ;
00354 
00355     switch (base_t) {
00356                 
00357         case T_COS :
00358         case T_SIN :
00359         case T_SIN_P :
00360         case T_SIN_I :
00361         case T_COS_I : 
00362         case T_COS_P : {
00363 
00364         for (int k=2 ; k<np+1 ; k++) {
00365             somt = 0 ;
00366             for (int j=0 ; j<nt ; j++) {
00367             somt += (*pi) * tab_theta(0, j, j0) ;
00368             pi++ ;  // theta suivant
00369             }
00370             resu += somt * tab_phi(k, k0) ;
00371         }
00372         break ;
00373         } 
00374         
00375         case T_COSSIN_C : 
00376         case T_COSSIN_S : 
00377         case T_COSSIN_SP : 
00378         case T_COSSIN_SI : 
00379         case T_COSSIN_CI : 
00380         case T_COSSIN_CP : {
00381 
00382         for (int k=2 ; k<np+1 ; k++) {
00383             int m_par = (k/2)%2 ;   // 0 pour m pair, 1 pour m impair
00384             somt = 0 ;
00385             for (int j=0 ; j<nt ; j++) {
00386             somt += (*pi) * tab_theta(m_par, j, j0) ;
00387             pi++ ;  // theta suivant
00388             }
00389             resu += somt * tab_phi(k, k0) ;
00390         }
00391         break ;
00392         }
00393      
00394         default: {
00395         cout << "Mtbl_cf::val_point_jk: unknown theta basis ! " << endl ;
00396         abort () ;
00397         }
00398            
00399     }   // fin des cas sur base_t 
00400     
00401     }   // fin du cas np > 1 
00402 
00403 
00404     // Menage
00405     delete [] coef_tp ;
00406     
00407     // Termine
00408     return resu ;  
00409  
00410 }
00411  
00412 
00413     //-------------------------------------------------------------//
00414     //      version for xi = 1                                 //
00415     //      and collocation point in (theta',phi')         //
00416     //-------------------------------------------------------------//
00417 
00418 double Mtbl_cf::val_out_bound_jk(int l, int j0, int k0) const {
00419 
00420 #ifndef NDEBUG
00421 // Bases de developpement : 
00422     int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
00423     assert((base_r == R_CHEB) || (base_r == R_CHEBU) || (base_r == R_CHEBP)
00424        || (base_r == R_CHEBI) || (base_r == R_CHEBPIM_P) || (base_r == R_CHEBPIM_I)
00425        || (base_r == R_CHEBPI_P) || (base_r == R_CHEBPI_I) || (base_r == R_JACO02)) ;
00426 #endif
00427 
00428     int nr = mg->get_nr(l) ;
00429     double resu = 0 ;
00430     for (int i=0; i<nr; i++)
00431     resu += operator()(l, k0, j0, i) ;
00432 
00433     return resu ;
00434 }
00435 
00436 
00437     //-------------------------------------------------------------//
00438     //      version for xi = -1                                //
00439     //      and collocation point in (theta',phi')         //
00440     //-------------------------------------------------------------//
00441 
00442 double Mtbl_cf::val_in_bound_jk(int l, int j0, int k0) const {
00443 
00444 #ifndef NDEBUG
00445 // Bases de developpement : 
00446     int base_r = (base.b[l] & MSQ_R) >> TRA_R ;
00447     assert((base_r == R_CHEB) || (base_r == R_CHEBU)) ;
00448 #endif
00449 
00450     int nr = mg->get_nr(l) ;
00451     double resu = 0 ;
00452     int pari = 1 ;
00453     for (int i=0; i<nr; i++) {
00454     resu += pari*operator()(l, k0, j0, i) ;
00455     pari = - pari ;
00456     }
00457 
00458     return resu ;
00459 }
00460