op_mult_ct.C

00001 /*
00002  *   Copyright (c) 1999-2001 Jerome Novak
00003  *
00004  *   This file is part of LORENE.
00005  *
00006  *   LORENE is free software; you can redistribute it and/or modify
00007  *   it under the terms of the GNU General Public License as published by
00008  *   the Free Software Foundation; either version 2 of the License, or
00009  *   (at your option) any later version.
00010  *
00011  *   LORENE is distributed in the hope that it will be useful,
00012  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  *   GNU General Public License for more details.
00015  *
00016  *   You should have received a copy of the GNU General Public License
00017  *   along with LORENE; if not, write to the Free Software
00018  *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00019  *
00020  */
00021 
00022 
00023 char op_mult_ct_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_ct.C,v 1.6 2009/10/09 14:00:54 j_novak Exp $" ;
00024 
00025 /*
00026  * $Id: op_mult_ct.C,v 1.6 2009/10/09 14:00:54 j_novak Exp $
00027  * $Log: op_mult_ct.C,v $
00028  * Revision 1.6  2009/10/09 14:00:54  j_novak
00029  * New bases T_cos and T_SIN.
00030  *
00031  * Revision 1.5  2007/12/21 10:43:37  j_novak
00032  * Corrected some bugs in the case nt=1 (1 point in theta).
00033  *
00034  * Revision 1.4  2007/10/05 12:37:20  j_novak
00035  * Corrected a few errors in the theta-nonsymmetric case (bases T_COSSIN_C and
00036  * T_COSSIN_S).
00037  *
00038  * Revision 1.3  2005/02/16 15:29:23  m_forot
00039  * Correct T_COSSIN_S and T_COSSIN_C cases
00040  *
00041  * Revision 1.2  2004/11/23 15:16:01  m_forot
00042  *
00043  * Added the bases for the cases without any equatorial symmetry
00044  *  (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
00045  *
00046  * Revision 1.1.1.1  2001/11/20 15:19:29  e_gourgoulhon
00047  * LORENE
00048  *
00049  * Revision 1.3  2000/02/24  16:41:51  eric
00050  * Initialisation a zero du tableau xo avant le calcul.
00051  *
00052  * Revision 1.2  1999/11/23  16:14:09  novak
00053  * *** empty log message ***
00054  *
00055  * Revision 1.1  1999/11/23  14:31:56  novak
00056  * Initial revision
00057  *
00058  *
00059  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_ct.C,v 1.6 2009/10/09 14:00:54 j_novak Exp $
00060  *
00061  */
00062 
00063 /* 
00064  * Ensemble des routines de base de multiplication par cost theta
00065  * (Utilisation interne)
00066  * 
00067  *  void _mult_ct_XXXX(Tbl * t, int & b)
00068  *  t   pointeur sur le Tbl d'entree-sortie
00069  *  b   base spectrale
00070  * 
00071  */
00072  
00073  
00074 // Fichier includes
00075 #include "tbl.h"
00076 
00077  
00078         //-----------------------------------
00079         // Routine pour les cas non prevus --
00080         //-----------------------------------
00081 
00082 void _mult_ct_pas_prevu(Tbl * tb, int& base) {
00083   cout << "mult_ct pas prevu..." << endl ;
00084   cout << "Tbl: " << tb << " base: " << base << endl ;
00085   abort () ;
00086   exit(-1) ;
00087 }
00088 
00089             //--------------
00090             // cas T_COS
00091             //--------------
00092             
00093 void _mult_ct_t_cos(Tbl* tb, int & b)
00094 {
00095 
00096   // Peut-etre rien a faire ?
00097   if (tb->get_etat() == ETATZERO) {
00098     int base_r = b & MSQ_R ;
00099     int base_p = b & MSQ_P ;
00100     switch(base_r){
00101     case(R_CHEBPI_P):
00102       b = R_CHEBPI_I | base_p | T_COS ;
00103       break ;
00104     case(R_CHEBPI_I):
00105       b = R_CHEBPI_P | base_p | T_COS ;
00106       break ;  
00107     default:
00108       b = base_r | base_p | T_COS ;
00109       break;
00110     }
00111     return ;
00112   }
00113     
00114   // Protection
00115   assert(tb->get_etat() == ETATQCQ) ;
00116   
00117   // Pour le confort : nbre de points reels.
00118   int nr = (tb->dim).dim[0] ;
00119   int nt = (tb->dim).dim[1] ;
00120   int np = (tb->dim).dim[2] ;
00121   np = np - 2 ;
00122   
00123   // zone de travail
00124   double* som = new double [nr] ;
00125   
00126   // pt. sur le tableau de double resultat
00127   double* xo = new double[(tb->dim).taille] ;
00128   
00129     // Initialisation a zero :
00130     for (int i=0; i<(tb->dim).taille; i++) {
00131     xo[i] = 0 ; 
00132     }
00133     
00134   // On y va...
00135   double* xi = tb->t ;
00136   double* xci = xi ;    // Pointeurs
00137   double* xco = xo ;    //  courants
00138   
00139   // k = 0
00140   
00141   // Dernier j: j = nt-1
00142   // Positionnement
00143   xci += nr * (nt-1) ;
00144   xco += nr * (nt-1) ;
00145   
00146   // Initialisation de som 
00147   for (int i=0 ; i<nr ; i++) {
00148     som[i] = 0.5*xci[i] ;
00149     xco[i] = 0. ;   // mise a zero du dernier coef en theta
00150   }
00151   
00152   // j suivants
00153   for (int j=nt-2 ; j > 0 ; j--) {
00154     // Positionnement
00155     xci -= 2*nr ;
00156     xco -= nr ;
00157     // Calcul
00158     for (int i=0 ; i<nr ; i++ ) {
00159       som[i] += 0.5*xci[i] ;
00160       xco[i] = som[i] ;
00161     }   // Fin de la boucle sur r
00162     xci += nr ;
00163     for (int i=0; i<nr; i++)
00164     som[i] = 0.5*xci[i] ;
00165   }   // Fin de la boucle sur theta
00166   // j = 0 : le facteur 2...
00167   xci -= nr ;
00168   for (int i=0 ; i<nr ; i++ ) {
00169       xco[i] += 0.5*xci[i] ;
00170   }
00171   xco -= nr ;
00172   for (int i=0; i<nr; i++)
00173       xco[i] = som[i] ;
00174 
00175   // Positionnement phi suivant
00176   xci += nr*nt ;
00177   xco += nr*nt ;
00178   
00179   // k = 1
00180   xci += nr*nt ;
00181   xco += nr*nt ;
00182   
00183   // k >= 2
00184   for (int k=2 ; k<np+1 ; k++) {
00185     
00186     // Dernier j: j = nt-1
00187     // Positionnement
00188     xci += nr * (nt-1) ;
00189     xco += nr * (nt-1) ;
00190     
00191     // Initialisation de som
00192     for (int i=0 ; i<nr ; i++) {
00193       som[i] = 0.5*xci[i] ;
00194       xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
00195     }
00196     
00197     // j suivants
00198     for (int j=nt-2 ; j > 0 ; j--) {
00199       // Positionnement
00200       xci -= 2*nr ;
00201       xco -= nr ;
00202       // Calcul
00203       for (int i=0 ; i<nr ; i++ ) {
00204     som[i] += 0.5*xci[i] ;
00205     xco[i] = som[i] ;
00206       } // Fin de la boucle sur r
00207       xci += nr ;
00208       for (int i=0; i<nr; i++)
00209       som[i] = 0.5*xci[i] ;
00210     }   // Fin de la boucle sur theta
00211     // j = 0 : le facteur 2...
00212     xci -= nr ;
00213     for (int i = 0; i<nr; i++) {
00214     xco[i] += 0.5*xci[i] ;
00215     }
00216     xco -= nr ;
00217     for (int i=0; i<nr; i++)
00218     xco[i] = som[i] ;
00219     // Positionnement phi suivant
00220     xci += nr*nt ;
00221     xco += nr*nt ;
00222   } // Fin de la boucle sur phi
00223   
00224   // On remet les choses la ou il faut
00225   delete [] tb->t ;
00226   tb->t = xo ;
00227   
00228   //Menage
00229   delete [] som ;
00230   
00231   // base de developpement
00232   int base_r = b & MSQ_R ;
00233   int base_p = b & MSQ_P ;
00234     switch(base_r){
00235     case(R_CHEBPI_P):
00236       b = R_CHEBPI_I | base_p | T_COS ;
00237       break ;
00238     case(R_CHEBPI_I):
00239       b = R_CHEBPI_P | base_p | T_COS ;
00240       break ;  
00241     default:
00242       b = base_r | base_p | T_COS ;
00243       break;
00244     }
00245 }
00246             
00247             //------------
00248             // cas T_SIN
00249             //------------
00250             
00251 void _mult_ct_t_sin(Tbl* tb, int & b)
00252 {
00253   // Peut-etre rien a faire ?
00254   if (tb->get_etat() == ETATZERO) {
00255     int base_r = b & MSQ_R ;
00256     int base_p = b & MSQ_P ;
00257     switch(base_r){
00258     case(R_CHEBPI_P):
00259       b = R_CHEBPI_I | base_p | T_SIN ;
00260       break ;
00261     case(R_CHEBPI_I):
00262       b = R_CHEBPI_P | base_p | T_SIN ;
00263       break ;  
00264     default:
00265       b = base_r | base_p | T_SIN ;
00266       break;
00267     }
00268     return ;
00269   }
00270   
00271   // Protection
00272   assert(tb->get_etat() == ETATQCQ) ;
00273   
00274   // Pour le confort : nbre de points reels.
00275   int nr = (tb->dim).dim[0] ;
00276   int nt = (tb->dim).dim[1] ;
00277   int np = (tb->dim).dim[2] ;
00278   np = np - 2 ;
00279   
00280   // zone de travail
00281   double* som = new double [nr] ;
00282   
00283   // pt. sur le tableau de double resultat
00284   double* xo = new double[(tb->dim).taille] ;
00285   
00286     // Initialisation a zero :
00287     for (int i=0; i<(tb->dim).taille; i++) {
00288     xo[i] = 0 ; 
00289     }
00290     
00291   // On y va...
00292   double* xi = tb->t ;
00293   double* xci = xi ;    // Pointeurs
00294   double* xco = xo ;    //  courants
00295   
00296   // k = 0
00297   
00298   // Dernier j: j = nt-1
00299   // Positionnement
00300   xci += nr * (nt-1) ;
00301   xco += nr * (nt-1) ;
00302   
00303   // Initialisation de som
00304   for (int i=0 ; i<nr ; i++) {
00305     som[i] = 0.5*xci[i] ;
00306     xco[i] = 0. ;
00307   }
00308   
00309   // j suivants
00310   for (int j=nt-2 ; j > 0 ; j--) {
00311     // Positionnement
00312     xci -= 2*nr ;
00313     xco -= nr ;
00314     // Calcul
00315     for (int i=0 ; i<nr ; i++ ) {
00316       som[i] += 0.5*xci[i] ;
00317       xco[i] = som[i] ;
00318     }   // Fin de la boucle sur r
00319     xci += nr ;
00320     for (int i=0; i<nr; i++) {
00321       som[i] = 0.5*xci[i] ;
00322     }
00323   }   // Fin de la boucle sur theta
00324   // j = 0 : sin(0*theta)
00325   xci -= nr ;
00326   xco -= nr  ;
00327   for (int i =0; i<nr ; i++) {
00328     xco[i] = 0.  ;
00329   }
00330   // Positionnement phi suivant
00331   xci += nr*nt ;
00332   xco += nr*nt ;
00333   
00334   // k = 1
00335   xci += nr*nt ;
00336   xco += nr*nt ;
00337   
00338   // k >= 2
00339   for (int k=2 ; k<np+1 ; k++) {
00340     
00341     // Dernier j: j = nt-1
00342     // Positionnement
00343     xci += nr * (nt-1) ;
00344     xco += nr * (nt-1) ;
00345     
00346     // Initialisation de som
00347     for (int i=0 ; i<nr ; i++) {
00348       som[i] = 0.5*xci[i] ;
00349       xco[i] = 0. ;
00350     }
00351     
00352     // j suivants
00353     for (int j=nt-2 ; j > 0 ; j--) {
00354       // Positionnement
00355       xci -= 2*nr ;
00356       xco -= nr ;
00357       // Calcul
00358       for (int i=0 ; i<nr ; i++ ) {
00359     som[i] += 0.5*xci[i] ;
00360     xco[i] = som[i] ;
00361       } // Fin de la boucle sur r
00362       xci += nr ;
00363       for (int i=0; i<nr; i++) {
00364     som[i] = 0.5*xci[i] ;
00365       }
00366     }   // Fin de la boucle sur theta
00367     // j = 0 : sin (0*theta) 
00368     xci -= nr ;
00369     xco -= nr ;
00370     for (int i=0; i<nr; i++) {
00371       xco[i] = 0.0 ;
00372     }
00373     // Positionnement phi suivant
00374     xci += nr*nt ;
00375     xco += nr*nt ;
00376   } // Fin de la boucle sur phi
00377   
00378   // On remet les choses la ou il faut
00379   delete [] tb->t ;
00380   tb->t = xo ;
00381   
00382   //Menage
00383   delete [] som ;
00384   
00385   // base de developpement
00386   int base_r = b & MSQ_R ;
00387   int base_p = b & MSQ_P ;
00388     switch(base_r){
00389     case(R_CHEBPI_P):
00390       b = R_CHEBPI_I | base_p | T_SIN ;
00391       break ;
00392     case(R_CHEBPI_I):
00393       b = R_CHEBPI_P | base_p | T_SIN ;
00394       break ;  
00395     default:
00396       b = base_r | base_p | T_SIN ;
00397       break;
00398     }  
00399 }
00400             //----------------
00401             // cas T_COS_P
00402             //----------------
00403             
00404 void _mult_ct_t_cos_p(Tbl* tb, int & b)
00405 {
00406 
00407   // Peut-etre rien a faire ?
00408   if (tb->get_etat() == ETATZERO) {
00409     int base_r = b & MSQ_R ;
00410     int base_p = b & MSQ_P ;
00411     b = base_r | base_p | T_COS_I ;
00412     return ;
00413   }
00414     
00415   // Protection
00416   assert(tb->get_etat() == ETATQCQ) ;
00417   
00418   // Pour le confort : nbre de points reels.
00419   int nr = (tb->dim).dim[0] ;
00420   int nt = (tb->dim).dim[1] ;
00421   int np = (tb->dim).dim[2] ;
00422   np = np - 2 ;
00423   
00424   // zone de travail
00425   double* som = new double [nr] ;
00426   
00427   // pt. sur le tableau de double resultat
00428   double* xo = new double[(tb->dim).taille] ;
00429   
00430     // Initialisation a zero :
00431     for (int i=0; i<(tb->dim).taille; i++) {
00432     xo[i] = 0 ; 
00433     }
00434     
00435   // On y va...
00436   double* xi = tb->t ;
00437   double* xci = xi ;    // Pointeurs
00438   double* xco = xo ;    //  courants
00439   
00440   // k = 0
00441   
00442   // Dernier j: j = nt-1
00443   // Positionnement
00444   xci += nr * (nt-1) ;
00445   xco += nr * (nt-1) ;
00446   
00447   // Initialisation de som 
00448   for (int i=0 ; i<nr ; i++) {
00449     som[i] = 0.5*xci[i] ;
00450     xco[i] = 0. ;   // mise a zero du dernier coef en theta
00451   }
00452   
00453   // j suivants
00454   for (int j=nt-2 ; j > 0 ; j--) {
00455     // Positionnement
00456     xci -= nr ;
00457     xco -= nr ;
00458     // Calcul
00459     for (int i=0 ; i<nr ; i++ ) {
00460       som[i] += 0.5*xci[i] ;
00461       xco[i] = som[i] ;
00462       som[i] = 0.5*xci[i] ;
00463     }   // Fin de la boucle sur r
00464   }   // Fin de la boucle sur theta
00465 
00466   if (nt > 1) {
00467       // j = 0 
00468       xci -= nr ;
00469       xco -= nr ;
00470       for (int i=0 ; i<nr ; i++ ) {
00471       xco[i] = som[i]+xci[i] ;
00472       } // Fin de la boucle sur r
00473   }
00474 
00475   // Positionnement phi suivant
00476   xci += nr*nt ;
00477   xco += nr*nt ;
00478   
00479   // k = 1
00480   xci += nr*nt ;
00481   xco += nr*nt ;
00482   
00483   // k >= 2
00484   for (int k=2 ; k<np+1 ; k++) {
00485     
00486     // Dernier j: j = nt-1
00487     // Positionnement
00488     xci += nr * (nt-1) ;
00489     xco += nr * (nt-1) ;
00490     
00491     // Initialisation de som
00492     for (int i=0 ; i<nr ; i++) {
00493       som[i] = 0.5*xci[i] ;
00494       xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
00495     }
00496     
00497     // j suivants
00498     for (int j=nt-2 ; j > 0 ; j--) {
00499       // Positionnement
00500       xci -= nr ;
00501       xco -= nr ;
00502       // Calcul
00503       for (int i=0 ; i<nr ; i++ ) {
00504     som[i] += 0.5*xci[i] ;
00505     xco[i] = som[i] ;
00506     som[i] = 0.5*xci[i] ;
00507       } // Fin de la boucle sur r
00508     }   // Fin de la boucle sur theta
00509     // j = 0
00510     xci -= nr ;
00511     xco -= nr ;
00512     for (int i = 0; i<nr; i++) {
00513       xco[i] = xci[i] + som[i] ;
00514     }
00515     // Positionnement phi suivant
00516     xci += nr*nt ;
00517     xco += nr*nt ;
00518   } // Fin de la boucle sur phi
00519   
00520   // On remet les choses la ou il faut
00521   delete [] tb->t ;
00522   tb->t = xo ;
00523   
00524   //Menage
00525   delete [] som ;
00526   
00527   // base de developpement
00528   int base_r = b & MSQ_R ;
00529   int base_p = b & MSQ_P ;
00530   b = base_r | base_p | T_COS_I ;
00531 }
00532             
00533             //--------------
00534             // cas T_SIN_P
00535             //--------------
00536             
00537 void _mult_ct_t_sin_p(Tbl* tb, int & b)
00538 {
00539   // Peut-etre rien a faire ?
00540   if (tb->get_etat() == ETATZERO) {
00541     int base_r = b & MSQ_R ;
00542     int base_p = b & MSQ_P ;
00543     b = base_r | base_p | T_SIN_I ;
00544     return ;
00545   }
00546   
00547   // Protection
00548   assert(tb->get_etat() == ETATQCQ) ;
00549   
00550   // Pour le confort : nbre de points reels.
00551   int nr = (tb->dim).dim[0] ;
00552   int nt = (tb->dim).dim[1] ;
00553   int np = (tb->dim).dim[2] ;
00554   np = np - 2 ;
00555   
00556   // zone de travail
00557   double* som = new double [nr] ;
00558   
00559   // pt. sur le tableau de double resultat
00560   double* xo = new double[(tb->dim).taille] ;
00561   
00562     // Initialisation a zero :
00563     for (int i=0; i<(tb->dim).taille; i++) {
00564     xo[i] = 0 ; 
00565     }
00566     
00567   // On y va...
00568   double* xi = tb->t ;
00569   double* xci = xi ;    // Pointeurs
00570   double* xco = xo ;    //  courants
00571   
00572   // k = 0
00573   
00574   // Dernier j: j = nt-1
00575   // Positionnement
00576   xci += nr * (nt-1) ;
00577   xco += nr * (nt-1) ;
00578   
00579   // Initialisation de som
00580   for (int i=0 ; i<nr ; i++) {
00581     som[i] = 0.5*xci[i] ;
00582     xco[i] = 0. ;
00583   }
00584   
00585   // j suivants
00586   for (int j=nt-2 ; j > 0 ; j--) {
00587     // Positionnement
00588     xci -= nr ;
00589     xco -= nr ;
00590     // Calcul
00591     for (int i=0 ; i<nr ; i++ ) {
00592       som[i] += 0.5*xci[i] ;
00593       xco[i] = som[i] ;
00594       som[i] = 0.5*xci[i] ;
00595     }   // Fin de la boucle sur r
00596   }   // Fin de la boucle sur theta
00597   // j = 0
00598   if (nt > 1) {
00599       xci -= nr ;
00600       xco -= nr  ;
00601       for (int i =0; i<nr ; i++) {
00602       xco[i] = som[i]  ;
00603       }
00604   }
00605   
00606   // Positionnement phi suivant
00607   xci += nr*nt ;
00608   xco += nr*nt ;
00609   
00610   // k = 1
00611   xci += nr*nt ;
00612   xco += nr*nt ;
00613   
00614   // k >= 2
00615   for (int k=2 ; k<np+1 ; k++) {
00616     
00617     // Dernier j: j = nt-1
00618     // Positionnement
00619     xci += nr * (nt-1) ;
00620     xco += nr * (nt-1) ;
00621     
00622     // Initialisation de som
00623     for (int i=0 ; i<nr ; i++) {
00624       som[i] = 0.5*xci[i] ;
00625       xco[i] = 0. ;
00626     }
00627     
00628     // j suivants
00629     for (int j=nt-2 ; j > 0 ; j--) {
00630       // Positionnement
00631       xci -= nr ;
00632       xco -= nr ;
00633       // Calcul
00634       for (int i=0 ; i<nr ; i++ ) {
00635     som[i] += 0.5*xci[i] ;
00636     xco[i] = som[i] ;
00637     som[i] = 0.5*xci[i] ;
00638       } // Fin de la boucle sur r
00639     }   // Fin de la boucle sur theta
00640     // j = 0 
00641     xci -= nr ;
00642     xco -= nr ;
00643     for (int i=0; i<nr; i++) {
00644       xco[i] = som[i] ;
00645     }
00646     // Positionnement phi suivant
00647     xci += nr*nt ;
00648     xco += nr*nt ;
00649   } // Fin de la boucle sur phi
00650   
00651   // On remet les choses la ou il faut
00652   delete [] tb->t ;
00653   tb->t = xo ;
00654   
00655   //Menage
00656   delete [] som ;
00657   
00658   // base de developpement
00659   int base_r = b & MSQ_R ;
00660   int base_p = b & MSQ_P ;
00661   b = base_r | base_p | T_SIN_I ;
00662   
00663 }
00664             //--------------
00665             // cas T_SIN_I
00666             //--------------
00667             
00668 void _mult_ct_t_sin_i(Tbl* tb, int & b)
00669 {
00670   // Peut-etre rien a faire ?
00671   if (tb->get_etat() == ETATZERO) {
00672     int base_r = b & MSQ_R ;
00673     int base_p = b & MSQ_P ;
00674     b = base_r | base_p | T_SIN_P ;
00675     return ;
00676   }
00677     
00678   // Protection
00679   assert(tb->get_etat() == ETATQCQ) ;
00680   
00681   // Pour le confort : nbre de points reels.
00682   int nr = (tb->dim).dim[0] ;
00683   int nt = (tb->dim).dim[1] ;
00684   int np = (tb->dim).dim[2] ;
00685   np = np - 2 ;
00686   
00687   // zone de travail
00688   double* som = new double [nr] ;
00689   
00690   // pt. sur le tableau de double resultat
00691   double* xo = new double[(tb->dim).taille] ;
00692   
00693     // Initialisation a zero :
00694     for (int i=0; i<(tb->dim).taille; i++) {
00695     xo[i] = 0 ; 
00696     }
00697     
00698   // On y va...
00699   double* xi = tb->t ;
00700   double* xci = xi ;    // Pointeurs
00701   double* xco = xo ;    //  courants
00702   
00703   // k = 0
00704   
00705   // Dernier j: j = nt-1
00706   // Positionnement
00707   xci += nr * (nt-1) ;
00708   xco += nr * (nt-1) ;
00709   
00710   // Initialisation de som 
00711   for (int i=0 ; i<nr ; i++) {
00712     som[i] = 0. ;
00713   }
00714   
00715   // j suivants
00716   for (int j=nt-1 ; j > 0 ; j--) {
00717     // Positionnement
00718     xci -= nr ;
00719     // Calcul
00720     for (int i=0 ; i<nr ; i++ ) {
00721       som[i] += 0.5*xci[i] ;
00722       xco[i] = som[i] ;
00723       som[i] = 0.5*xci[i] ;
00724     }   // Fin de la boucle sur r
00725     xco -= nr ;
00726   }   // Fin de la boucle sur theta
00727   for (int i=0; i<nr; i++) {
00728     xco[i] = 0. ;
00729   }
00730 
00731   // Positionnement phi suivant
00732   xci += nr*nt ;
00733   xco += nr*nt ;
00734   
00735   // k = 1
00736   xci += nr*nt ;
00737   xco += nr*nt ;
00738     
00739   // k >= 2
00740   for (int k=2 ; k<np+1 ; k++) {
00741     
00742     // Dernier j: j = nt-1
00743     // Positionnement
00744     xci += nr * (nt-1) ;
00745     xco += nr * (nt-1) ;
00746     
00747     // Initialisation de som
00748     for (int i=0 ; i<nr ; i++) {
00749       som[i] = 0. ;
00750     }
00751     
00752     // j suivants
00753     for (int j=nt-1 ; j > 0 ; j--) {
00754       // Positionnement
00755       xci -= nr ;
00756       // Calcul
00757       for (int i=0 ; i<nr ; i++ ) {
00758     som[i] += 0.5*xci[i] ;
00759     xco[i] = som[i] ;
00760     som[i] = 0.5*xci[i] ;
00761       } // Fin de la boucle sur r
00762       xco -= nr ;
00763     }   // Fin de la boucle sur theta
00764   for (int i=0; i<nr; i++) {
00765     xco[i] = 0. ;
00766   }
00767 
00768     // Positionnement phi suivant
00769     xci += nr*nt ;
00770     xco += nr*nt ;
00771   } // Fin de la boucle sur phi
00772 
00773   // On remet les choses la ou il faut
00774   delete [] tb->t ;
00775   tb->t = xo ;
00776   
00777   //Menage
00778   delete [] som ;
00779   
00780   // base de developpement
00781   int base_r = b & MSQ_R ;
00782   int base_p = b & MSQ_P ;
00783   b = base_r | base_p | T_SIN_P ;
00784   
00785 }
00786                     //---------------------
00787                     // cas T_COS_I
00788             //---------------------
00789 void _mult_ct_t_cos_i(Tbl* tb, int & b)
00790 {
00791   // Peut-etre rien a faire ?
00792   if (tb->get_etat() == ETATZERO) {
00793     int base_r = b & MSQ_R ;
00794     int base_p = b & MSQ_P ;
00795     b = base_r | base_p | T_COS_P ;
00796     return ;
00797   }
00798   
00799   // Protection
00800   assert(tb->get_etat() == ETATQCQ) ;
00801   
00802   // Pour le confort : nbre de points reels.
00803   int nr = (tb->dim).dim[0] ;
00804   int nt = (tb->dim).dim[1] ;
00805   int np = (tb->dim).dim[2] ;
00806   np = np - 2 ;
00807   
00808   // zone de travail
00809   double* som = new double [nr] ;
00810   
00811   // pt. sur le tableau de double resultat
00812   double* xo = new double[(tb->dim).taille] ;
00813   
00814     // Initialisation a zero :
00815     for (int i=0; i<(tb->dim).taille; i++) {
00816     xo[i] = 0 ; 
00817     }
00818     
00819   // On y va...
00820   double* xi = tb->t ;
00821   double* xci = xi ;    // Pointeurs
00822   double* xco = xo ;    //  courants
00823   
00824   // k = 0
00825   
00826   // Dernier j: j = nt-1
00827   // Positionnement
00828   xci += nr * (nt-1) ;
00829   xco += nr * (nt-1) ;
00830   
00831   // Initialisation de som 
00832   for (int i=0 ; i<nr ; i++) {
00833     som[i] = 0. ;
00834   }
00835   
00836   // j suivants
00837   for (int j=nt-1 ; j > 0 ; j--) {
00838     // Positionnement
00839     xci -= nr ;
00840     // Calcul
00841     for (int i=0 ; i<nr ; i++ ) {
00842       som[i] += 0.5*xci[i] ;
00843       xco[i] = som[i] ;
00844       som[i] = 0.5*xci[i] ;
00845     }   // Fin de la boucle sur r
00846     xco -= nr ;
00847   }   // Fin de la boucle sur theta
00848       // premier theta
00849   for (int i=0 ; i<nr ; i++) {
00850       xco[i] = som[i] ;
00851   }
00852   // Positionnement phi suivant
00853   xci += nr*nt ;
00854   xco += nr*nt ;
00855   
00856   // k = 1
00857   xci += nr*nt ;
00858   xco += nr*nt ;
00859     
00860   // k >= 2
00861   for (int k=2 ; k<np+1 ; k++) {
00862     
00863     // Dernier j: j = nt-1
00864     // Positionnement
00865     xci += nr * (nt-1) ;
00866     xco += nr * (nt-1) ;
00867     
00868     // Initialisation de som
00869     for (int i=0 ; i<nr ; i++) {
00870       som[i] = 0. ;
00871     }
00872     
00873     // j suivants
00874     for (int j=nt-1 ; j > 0 ; j--) {
00875       // Positionnement
00876       xci -= nr ;
00877       // Calcul
00878       for (int i=0 ; i<nr ; i++ ) {
00879     som[i] += 0.5*xci[i] ;
00880     xco[i] = som[i] ;
00881     som[i] = 0.5*xci[i] ;
00882       } // Fin de la boucle sur r
00883       xco -= nr ;
00884     }   // Fin de la boucle sur theta
00885     // premier theta
00886     for (int i=0 ; i<nr ; i++) {
00887       xco[i] = som[i] ;
00888     }
00889     // Positionnement phi suivant
00890     xci += nr*nt ;
00891     xco += nr*nt ;
00892   } // Fin de la boucle sur phi
00893   
00894   // On remet les choses la ou il faut
00895   delete [] tb->t ;
00896   tb->t = xo ;
00897   
00898   //Menage
00899   delete [] som ;
00900   
00901   // base de developpement
00902   int base_r = b & MSQ_R ;
00903   int base_p = b & MSQ_P ;
00904   b = base_r | base_p | T_COS_P ;
00905   
00906 }
00907             //----------------------
00908                         // cas T_COSSIN_CP
00909             //----------------------
00910 void _mult_ct_t_cossin_cp(Tbl* tb, int & b)
00911 {
00912   // Peut-etre rien a faire ?
00913   if (tb->get_etat() == ETATZERO) {
00914     int base_r = b & MSQ_R ;
00915     int base_p = b & MSQ_P ;
00916     b = base_r | base_p | T_COSSIN_CI ;
00917     return ;
00918   }
00919   
00920   // Protection
00921   assert(tb->get_etat() == ETATQCQ) ;
00922   
00923   // Pour le confort : nbre de points reels.
00924   int nr = (tb->dim).dim[0] ;
00925   int nt = (tb->dim).dim[1] ;
00926   int np = (tb->dim).dim[2] ;
00927   np = np - 2 ;
00928     
00929   // zone de travail
00930   double* som = new double [nr] ;
00931   
00932   // pt. sur le tableau de double resultat
00933   double* xo = new double[(tb->dim).taille] ;
00934     
00935     // Initialisation a zero :
00936     for (int i=0; i<(tb->dim).taille; i++) {
00937     xo[i] = 0 ; 
00938     }
00939     
00940   // On y va...
00941   double* xi = tb->t ;
00942   double* xci = xi ;    // Pointeurs
00943   double* xco = xo ;    //  courants
00944   
00945   // k = 0
00946   int m = 0 ;       // Parite de l'harmonique en phi
00947   // Dernier j: j = nt-1
00948   // Positionnement
00949   xci += nr * (nt-1) ;
00950   xco += nr * (nt-1) ;
00951     
00952   // Initialisation de som
00953   for (int i=0 ; i<nr ; i++) {
00954     som[i] = 0.5*xci[i] ;
00955     xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
00956   }
00957     
00958   // j suivants
00959   for (int j=nt-2 ; j > 0 ; j--) {
00960     // Positionnement
00961     xci -= nr ;
00962     xco -= nr ;
00963     // Calcul
00964     for (int i=0 ; i<nr ; i++ ) {
00965       som[i] += 0.5*xci[i] ;
00966       xco[i] = som[i] ;
00967       som[i] = 0.5*xci[i] ;
00968     }   // Fin de la boucle sur r
00969   }   // Fin de la boucle sur theta
00970   // j = 0
00971   xci -= nr ;
00972   xco -= nr ;
00973   for (int i = 0; i<nr; i++) {
00974     xco[i] = xci[i] + som[i] ;
00975   }
00976   // Positionnement phi suivant
00977   xci += nr*nt ;
00978   xco += nr*nt ;
00979   
00980   // k=1
00981   xci += nr*nt ;
00982   xco += nr*nt ;
00983   
00984   // k >= 2
00985   for (int k=2 ; k<np+1 ; k++) {
00986     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
00987     
00988     switch(m) {
00989     case 0:     // m pair: cos(pair)
00990       // Dernier j: j = nt-1
00991       // Positionnement
00992       xci += nr * (nt-1) ;
00993       xco += nr * (nt-1) ;
00994       
00995       // Initialisation de som
00996       for (int i=0 ; i<nr ; i++) {
00997     som[i] = 0.5*xci[i] ;
00998     xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
00999       }
01000       
01001       // j suivants
01002       for (int j=nt-2 ; j > 0 ; j--) {
01003     // Positionnement
01004     xci -= nr ;
01005     xco -= nr ;
01006     // Calcul
01007     for (int i=0 ; i<nr ; i++ ) {
01008       som[i] += 0.5*xci[i] ;
01009       xco[i] = som[i] ;
01010       som[i] = 0.5*xci[i] ;
01011     }   // Fin de la boucle sur r
01012       }   // Fin de la boucle sur theta
01013       // j = 0
01014       xci -= nr ;
01015       xco -= nr ;
01016       for (int i = 0; i<nr; i++) {
01017     xco[i] = xci[i] + som[i] ;
01018       }
01019       // Positionnement phi suivant
01020       xci += nr*nt ;
01021       xco += nr*nt ;
01022       break ;
01023       
01024     case 1:     // m impair: sin(impair)
01025       // Dernier j: j = nt-1
01026       // Positionnement
01027       xci += nr * (nt-1) ;
01028       xco += nr * (nt-1) ;
01029       
01030       // Initialisation de som
01031       for (int i=0 ; i<nr ; i++) {
01032     som[i] = 0. ;
01033       }
01034     
01035       // j suivants
01036       for (int j=nt-1 ; j > 0 ; j--) {
01037     // Positionnement
01038     xci -= nr ;
01039     // Calcul
01040     for (int i=0 ; i<nr ; i++ ) {
01041       som[i] += 0.5*xci[i] ;
01042       xco[i] = som[i] ;
01043       som[i] = 0.5*xci[i] ;
01044     }   // Fin de la boucle sur r
01045     xco -= nr ;
01046       }   // Fin de la boucle sur theta
01047       for (int i=0; i<nr; i++) {
01048     xco[i] = 0. ;
01049       }
01050 
01051       // Positionnement phi suivant
01052       xci += nr*nt ;
01053       xco += nr*nt ;
01054       break;
01055     } // Fin du switch
01056   } // Fin de la boucle sur phi
01057   
01058   // On remet les choses la ou il faut
01059   delete [] tb->t ;
01060   tb->t = xo ;
01061   
01062   //Menage
01063   delete [] som ;
01064   
01065   // base de developpement
01066   int base_r = b & MSQ_R ;
01067   int base_p = b & MSQ_P ;
01068   b = base_r | base_p | T_COSSIN_CI ;
01069 }
01070             
01071             //------------------
01072             // cas T_COSSIN_CI
01073             //----------------
01074 void _mult_ct_t_cossin_ci(Tbl* tb, int & b)
01075 {
01076   // Peut-etre rien a faire ?
01077   if (tb->get_etat() == ETATZERO) {
01078     int base_r = b & MSQ_R ;
01079     int base_p = b & MSQ_P ;
01080     b = base_r | base_p | T_COSSIN_CP ;
01081     return ;
01082   }
01083   
01084   // Protection
01085   assert(tb->get_etat() == ETATQCQ) ;
01086   
01087   // Pour le confort : nbre de points reels.
01088   int nr = (tb->dim).dim[0] ;
01089   int nt = (tb->dim).dim[1] ;
01090   int np = (tb->dim).dim[2] ;
01091   np = np - 2 ;
01092   
01093   // zone de travail
01094   double* som = new double [nr] ;
01095   
01096   // pt. sur le tableau de double resultat
01097   double* xo = new double[(tb->dim).taille] ;
01098   
01099     // Initialisation a zero :
01100     for (int i=0; i<(tb->dim).taille; i++) {
01101     xo[i] = 0 ; 
01102     }
01103     
01104   // On y va...
01105   double* xi = tb->t ;
01106   double* xci = xi ;    // Pointeurs
01107   double* xco = xo ;    //  courants
01108   
01109   // k = 0
01110   int m = 0 ;       // Parite de l'harmonique en phi
01111   // Dernier j: j = nt-1
01112   // Positionnement
01113   xci += nr * (nt-1) ;
01114   xco += nr * (nt-1) ;
01115   
01116   // Initialisation de som
01117   for (int i=0 ; i<nr ; i++) {
01118     som[i] = 0. ;
01119   }
01120   
01121   // j suivants
01122   for (int j=nt-1 ; j > 0 ; j--) {
01123     // Positionnement
01124     xci -= nr ;
01125     // Calcul
01126     for (int i=0 ; i<nr ; i++ ) {
01127       som[i] += 0.5*xci[i] ;
01128       xco[i] = som[i] ;
01129       som[i] = 0.5*xci[i] ;
01130     }   // Fin de la boucle sur r
01131     xco -= nr ;
01132   }   // Fin de la boucle sur theta
01133   // premier theta
01134   for (int i=0 ; i<nr ; i++) {
01135     xco[i] = som[i] ;
01136   }
01137   // Positionnement phi suivant
01138   xci += nr*nt ;
01139   xco += nr*nt ;
01140   
01141   // k=1
01142   xci += nr*nt ;
01143   xco += nr*nt ;
01144   
01145   // k >= 2
01146   for (int k=2 ; k<np+1 ; k++) {
01147     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
01148     
01149     switch(m) {
01150     case 0:     // m pair: cos(impair)
01151       // Dernier j: j = nt-1
01152       // Positionnement
01153       xci += nr * (nt-1) ;
01154       xco += nr * (nt-1) ;
01155       
01156       // Initialisation de som
01157       for (int i=0 ; i<nr ; i++) {
01158     som[i] = 0. ;
01159       }
01160       
01161       // j suivants
01162       for (int j=nt-1 ; j > 0 ; j--) {
01163     // Positionnement
01164     xci -= nr ;
01165     // Calcul
01166     for (int i=0 ; i<nr ; i++ ) {
01167       som[i] += 0.5*xci[i] ;
01168       xco[i] = som[i] ;
01169       som[i] = 0.5*xci[i] ;
01170     }   // Fin de la boucle sur r
01171     xco -= nr ;
01172       }   // Fin de la boucle sur theta
01173       // premier theta
01174       for (int i=0 ; i<nr ; i++) {
01175     xco[i] = som[i] ;
01176       }
01177       // Positionnement phi suivant
01178       xci += nr*nt ;
01179       xco += nr*nt ;
01180       break ;
01181       
01182     case 1:
01183       // Dernier j: j = nt-1
01184       // Positionnement
01185       xci += nr * (nt-1) ;
01186       xco += nr * (nt-1) ;
01187       
01188       // Initialisation de som
01189       for (int i=0 ; i<nr ; i++) {
01190     som[i] = 0.5*xci[i] ;
01191     xco[i] = 0. ;
01192       }
01193       
01194       // j suivants
01195       for (int j=nt-2 ; j > 0 ; j--) {
01196     // Positionnement
01197     xci -= nr ;
01198     xco -= nr ;
01199     // Calcul
01200     for (int i=0 ; i<nr ; i++ ) {
01201       som[i] += 0.5*xci[i] ;
01202       xco[i] = som[i] ;
01203       som[i] = 0.5*xci[i] ;
01204     }   // Fin de la boucle sur r
01205       }   // Fin de la boucle sur theta
01206       // j = 0 
01207       xci -= nr ;
01208       xco -= nr ;
01209       for (int i=0; i<nr; i++) {
01210     xco[i] = som[i] ;
01211       }
01212       // Positionnement phi suivant
01213       xci += nr*nt ;
01214       xco += nr*nt ;
01215       break ;
01216     }   // Fin du switch
01217   }   // Fin de la boucle en phi
01218   
01219   // On remet les choses la ou il faut
01220   delete [] tb->t ;
01221   tb->t = xo ;
01222   
01223   //Menage
01224   delete [] som ;
01225   
01226   // base de developpement
01227   int base_r = b & MSQ_R ;
01228   int base_p = b & MSQ_P ;
01229   b = base_r | base_p | T_COSSIN_CP ;
01230 }
01231 
01232             //--------------------- 
01233             // cas T_COSSIN_SI
01234             //----------------
01235 void _mult_ct_t_cossin_si(Tbl* tb, int & b)
01236 {
01237   // Peut-etre rien a faire ?
01238   if (tb->get_etat() == ETATZERO) {
01239     int base_r = b & MSQ_R ;
01240     int base_p = b & MSQ_P ;
01241     b = base_r | base_p | T_COSSIN_SP ;
01242     return ;
01243   }
01244   
01245   // Protection
01246   assert(tb->get_etat() == ETATQCQ) ;
01247   
01248   // Pour le confort : nbre de points reels.
01249   int nr = (tb->dim).dim[0] ;
01250   int nt = (tb->dim).dim[1] ;
01251   int np = (tb->dim).dim[2] ;
01252   np = np - 2 ;
01253   
01254   // zone de travail
01255   double* som = new double [nr] ;
01256   
01257   // pt. sur le tableau de double resultat
01258   double* xo = new double[(tb->dim).taille] ;
01259   
01260     // Initialisation a zero :
01261     for (int i=0; i<(tb->dim).taille; i++) {
01262     xo[i] = 0 ; 
01263     }
01264     
01265   // On y va...
01266   double* xi = tb->t ;
01267   double* xci = xi ;    // Pointeurs
01268   double* xco = xo ;    //  courants
01269   
01270   // k = 0
01271   int m = 0 ;       // Parite de l'harmonique en phi
01272   // Dernier j: j = nt-1
01273   // Positionnement
01274   xci += nr * (nt-1) ;
01275   xco += nr * (nt-1) ;
01276   
01277   // Initialisation de som
01278   for (int i=0 ; i<nr ; i++) {
01279     som[i] = 0. ;
01280   }
01281   
01282   // j suivants
01283   for (int j=nt-1 ; j > 0 ; j--) {
01284     // Positionnement
01285     xci -= nr ;
01286     // Calcul
01287     for (int i=0 ; i<nr ; i++ ) {
01288       som[i] += 0.5*xci[i] ;
01289       xco[i] = som[i] ;
01290       som[i] = 0.5*xci[i] ;
01291     }   // Fin de la boucle sur r
01292     xco -= nr ;
01293   }   // Fin de la boucle sur theta
01294   for (int i=0; i<nr; i++) {
01295     xco[i] = 0. ;
01296   }
01297 
01298   // Positionnement phi suivant
01299   xci += nr*nt ;
01300   xco += nr*nt ;
01301   
01302   // k=1 
01303   xci += nr*nt ;
01304   xco += nr*nt ;
01305   
01306   // k >= 2
01307   for (int k=2 ; k<np+1 ; k++) {
01308     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
01309     
01310     switch(m) {
01311     case 0:     // m pair: sin(impair)
01312       // Dernier j: j = nt-1
01313       // Positionnement
01314       xci += nr * (nt-1) ;
01315       xco += nr * (nt-1) ;
01316     
01317       // Initialisation de som
01318       for (int i=0 ; i<nr ; i++) {
01319     som[i] = 0. ;
01320       }
01321       
01322       // j suivants
01323       for (int j=nt-1 ; j > 0 ; j--) {
01324     // Positionnement
01325     xci -= nr ;
01326         // Calcul
01327     for (int i=0 ; i<nr ; i++ ) {
01328       som[i] += 0.5*xci[i] ;
01329       xco[i] = som[i] ;
01330       som[i] = 0.5*xci[i] ;
01331     }   // Fin de la boucle sur r
01332     xco -= nr ;
01333       }   // Fin de la boucle sur theta
01334       for (int i=0; i<nr; i++) {
01335     xco[i] = 0. ;
01336       }
01337 
01338       // Positionnement phi suivant
01339       xci += nr*nt ;
01340       xco += nr*nt ;
01341       break ;
01342       
01343     case 1: // m impair cos(pair)
01344       // Dernier j: j = nt-1
01345       // Positionnement
01346       xci += nr * (nt-1) ;
01347       xco += nr * (nt-1) ;
01348       
01349       // Initialisation de som
01350       for (int i=0 ; i<nr ; i++) {
01351     som[i] = 0.5*xci[i] ;
01352     xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
01353       }
01354       
01355       // j suivants
01356       for (int j=nt-2 ; j > 0 ; j--) {
01357     // Positionnement
01358     xci -= nr ;
01359     xco -= nr ;
01360     // Calcul
01361     for (int i=0 ; i<nr ; i++ ) {
01362       som[i] += 0.5*xci[i] ;
01363       xco[i] = som[i] ;
01364       som[i] = 0.5*xci[i] ;
01365     }   // Fin de la boucle sur r
01366       }   // Fin de la boucle sur theta
01367       // j = 0
01368       xci -= nr ;
01369       xco -= nr ;
01370       for (int i = 0; i<nr; i++) {
01371     xco[i] = xci[i] + som[i] ;
01372       }
01373       // Positionnement phi suivant
01374       xci += nr*nt ;
01375       xco += nr*nt ;
01376       break ;
01377     }   // Fin du switch
01378   }   // Fin de la boucle en phi
01379   
01380   // On remet les choses la ou il faut
01381   delete [] tb->t ;
01382   tb->t = xo ;
01383   
01384   //Menage
01385   delete [] som ;
01386   
01387   // base de developpement
01388   int base_r = b & MSQ_R ;
01389   int base_p = b & MSQ_P ;
01390   b = base_r | base_p | T_COSSIN_SP ;
01391   
01392   
01393 }
01394             //--------------------- 
01395             // cas T_COSSIN_SP
01396             //----------------
01397 void _mult_ct_t_cossin_sp(Tbl* tb, int & b)
01398 {
01399   // Peut-etre rien a faire ?
01400   if (tb->get_etat() == ETATZERO) {
01401     int base_r = b & MSQ_R ;
01402     int base_p = b & MSQ_P ;
01403     b = base_r | base_p | T_COSSIN_SI ;
01404     return ;
01405   }
01406   
01407   // Protection
01408   assert(tb->get_etat() == ETATQCQ) ;
01409   
01410   // Pour le confort : nbre de points reels.
01411   int nr = (tb->dim).dim[0] ;
01412   int nt = (tb->dim).dim[1] ;
01413   int np = (tb->dim).dim[2] ;
01414   np = np - 2 ;
01415   
01416   // zone de travail
01417   double* som = new double [nr] ;
01418   
01419   // pt. sur le tableau de double resultat
01420   double* xo = new double[(tb->dim).taille] ;
01421   
01422     // Initialisation a zero :
01423     for (int i=0; i<(tb->dim).taille; i++) {
01424     xo[i] = 0 ; 
01425     }
01426     
01427   // On y va...
01428   double* xi = tb->t ;
01429   double* xci = xi ;    // Pointeurs
01430   double* xco = xo ;    //  courants
01431   
01432   // k = 0
01433   int m = 0 ;       // Parite de l'harmonique en phi
01434   // Dernier j: j = nt-1
01435   // Positionnement
01436   xci += nr * (nt-1) ;
01437   xco += nr * (nt-1) ;
01438     
01439   // Initialisation de som
01440   for (int i=0 ; i<nr ; i++) {
01441     som[i] = 0.5*xci[i] ;
01442     xco[i] = 0. ;
01443   }
01444     
01445   // j suivants
01446   for (int j=nt-2 ; j > 0 ; j--) {
01447     // Positionnement
01448     xci -= nr ;
01449     xco -= nr ;
01450     // Calcul
01451     for (int i=0 ; i<nr ; i++ ) {
01452       som[i] += 0.5*xci[i] ;
01453       xco[i] = som[i] ;
01454       som[i] = 0.5*xci[i] ;
01455     }   // Fin de la boucle sur r
01456   }   // Fin de la boucle sur theta
01457   // j = 0 
01458   xci -= nr ;
01459   xco -= nr ;
01460   for (int i=0; i<nr; i++) {
01461     xco[i] = som[i] ;
01462   }
01463   // Positionnement phi suivant
01464   xci += nr*nt ;
01465   xco += nr*nt ;
01466   
01467   // k=1 
01468   xci += nr*nt ;
01469   xco += nr*nt ;
01470   
01471   for (int k=2 ; k<np+1 ; k++) {
01472     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
01473     
01474     switch(m) {
01475     case 1:     // m impair: cos(impair)
01476       // Dernier j: j = nt-1
01477       // Positionnement
01478       xci += nr * (nt-1) ;
01479       xco += nr * (nt-1) ;
01480       
01481       // Initialisation de som
01482       for (int i=0 ; i<nr ; i++) {
01483     som[i] = 0. ;
01484       }
01485       
01486       // j suivants
01487       for (int j=nt-1 ; j > 0 ; j--) {
01488     // Positionnement
01489     xci -= nr ;
01490     // Calcul
01491     for (int i=0 ; i<nr ; i++ ) {
01492       som[i] += 0.5*xci[i] ;
01493       xco[i] = som[i] ;
01494       som[i] = 0.5*xci[i] ;
01495     }   // Fin de la boucle sur r
01496     xco -= nr ;
01497       }   // Fin de la boucle sur theta
01498       // premier theta
01499       for (int i=0 ; i<nr ; i++) {
01500     xco[i] = som[i] ;
01501       }
01502       // Positionnement phi suivant
01503       xci += nr*nt ;
01504       xco += nr*nt ;
01505       break ;
01506       
01507     case 0:     // m pair: sin(pair)
01508       // Dernier j: j = nt-1
01509       // Positionnement
01510       xci += nr * (nt-1) ;
01511       xco += nr * (nt-1) ;
01512       
01513       // Initialisation de som
01514       for (int i=0 ; i<nr ; i++) {
01515     som[i] = 0.5*xci[i] ;
01516     xco[i] = 0. ;
01517       }
01518       
01519       // j suivants
01520       for (int j=nt-2 ; j > 0 ; j--) {
01521     // Positionnement
01522     xci -= nr ;
01523     xco -= nr ;
01524     // Calcul
01525     for (int i=0 ; i<nr ; i++ ) {
01526       som[i] += 0.5*xci[i] ;
01527       xco[i] = som[i] ;
01528       som[i] = 0.5*xci[i] ;
01529     }   // Fin de la boucle sur r
01530       }   // Fin de la boucle sur theta
01531       // j = 0 
01532       xci -= nr ;
01533       xco -= nr ;
01534       for (int i=0; i<nr; i++) {
01535     xco[i] = som[i] ;
01536       }
01537       // Positionnement phi suivant
01538       xci += nr*nt ;
01539       xco += nr*nt ;
01540       break ;
01541     }   // Fin du switch
01542   }   // Fin de la boucle en phi
01543   
01544   // On remet les choses la ou il faut
01545   delete [] tb->t ;
01546   tb->t = xo ;
01547   
01548   //Menage
01549   delete [] som ;
01550   
01551   // base de developpement
01552   int base_r = b & MSQ_R ;
01553   int base_p = b & MSQ_P ;
01554   b = base_r | base_p | T_COSSIN_SI ;
01555     
01556 }
01557 
01558             //----------------------
01559                         // cas T_COSSIN_C
01560             //----------------------
01561 void _mult_ct_t_cossin_c(Tbl* tb, int & b)
01562 {
01563   // Peut-etre rien a faire ?
01564   if (tb->get_etat() == ETATZERO) {
01565     int base_r = b & MSQ_R ;
01566     int base_p = b & MSQ_P ;
01567     switch(base_r){
01568     case(R_CHEBPI_P):
01569       b = R_CHEBPI_I | base_p | T_COSSIN_C ;
01570       break ;
01571     case(R_CHEBPI_I):
01572       b = R_CHEBPI_P | base_p | T_COSSIN_C ;
01573       break ;  
01574     default:
01575       b = base_r | base_p | T_COSSIN_C ;
01576       break;
01577     }
01578     return ;
01579   }
01580   
01581   // Protection
01582   assert(tb->get_etat() == ETATQCQ) ;
01583   
01584   // Pour le confort : nbre de points reels.
01585   int nr = (tb->dim).dim[0] ;
01586   int nt = (tb->dim).dim[1] ;
01587   int np = (tb->dim).dim[2] ;
01588   np = np - 2 ;
01589     
01590   // zone de travail
01591   double* som = new double [nr] ;
01592   
01593   // pt. sur le tableau de double resultat
01594   double* xo = new double[(tb->dim).taille] ;
01595     
01596     // Initialisation a zero :
01597     for (int i=0; i<(tb->dim).taille; i++) {
01598     xo[i] = 0 ; 
01599     }
01600     
01601   // On y va...
01602   double* xi = tb->t ;
01603   double* xci = xi ;    // Pointeurs
01604   double* xco = xo ;    //  courants
01605   
01606   // k = 0
01607   int m = 0 ;       // Parite de l'harmonique en phi
01608   // Dernier j: j = nt-1
01609   // Positionnement
01610   xci += nr * (nt-1) ;
01611   xco += nr * (nt-1) ;
01612     
01613   // Initialisation de som
01614   for (int i=0 ; i<nr ; i++) {
01615     som[i] = 0.5*xci[i] ;
01616     xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
01617   }
01618     
01619   // j suivants
01620   for (int j=nt-2 ; j > 0 ; j--) {
01621     // Positionnement
01622     xci -= 2*nr ;
01623     xco -= nr ;
01624     // Calcul
01625     for (int i=0 ; i<nr ; i++ ) {
01626       som[i] += 0.5*xci[i] ;
01627       xco[i] = som[i] ;
01628     }   // Fin de la boucle sur r
01629     xci += nr ;
01630     for (int i=0 ; i<nr ; i++ ) {
01631       som[i] = 0.5*xci[i] ;
01632     }
01633   }   // Fin de la boucle sur theta
01634   // j = 0 : le facteur 2...
01635   xci -= nr ;
01636   for (int i=0; i<nr; i++) {
01637       xco[i] += 0.5*xci[i] ;
01638   }
01639   xco -= nr ;
01640   for (int i = 0; i<nr; i++) {
01641     xco[i] = som[i] ;
01642   }
01643   // Positionnement phi suivant
01644   xci += nr*nt ;
01645   xco += nr*nt ;
01646   
01647   // k=1
01648   xci += nr*nt ;
01649   xco += nr*nt ;
01650   
01651   // k >= 2
01652   for (int k=2 ; k<np+1 ; k++) {
01653     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
01654     
01655     switch(m) {
01656     case 0:     // m pair: cos
01657       // Dernier j: j = nt-1
01658       // Positionnement
01659       xci += nr * (nt-1) ;
01660       xco += nr * (nt-1) ;
01661       
01662       // Initialisation de som
01663       for (int i=0 ; i<nr ; i++) {
01664     som[i] = 0.5*xci[i] ;
01665     xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
01666       }
01667       
01668       // j suivants
01669       for (int j=nt-2 ; j > 0 ; j--) {
01670     // Positionnement
01671     xci -= 2*nr ;
01672     xco -= nr ;
01673     // Calcul
01674     for (int i=0 ; i<nr ; i++ ) {
01675       som[i] += 0.5*xci[i] ;
01676       xco[i] = som[i] ;
01677     }   // Fin de la boucle sur r
01678     xci += nr ;
01679     for (int i=0 ; i<nr ; i++ ) {
01680       som[i] = 0.5*xci[i] ;
01681     }
01682       }   // Fin de la boucle sur theta
01683       // j = 0 : le facteur 2...
01684       xci -= nr ;
01685       for (int i=0; i<nr; i++) {
01686       xco[i] += 0.5*xci[i] ;
01687       }
01688       xco -= nr ;
01689       for (int i = 0; i<nr; i++) {
01690     xco[i] = som[i] ;
01691       }
01692       // Positionnement phi suivant
01693       xci += nr*nt ;
01694       xco += nr*nt ;
01695       break ;
01696       
01697     case 1:     // m impair: sin
01698       // Dernier j: j = nt-1
01699       // Positionnement
01700       xci += nr * (nt-1) ;
01701       xco += nr * (nt-1) ;
01702       
01703       // Initialisation de som
01704       for (int i=0 ; i<nr ; i++) {
01705     som[i] = 0.5*xci[i] ;
01706     xco[i] = 0.0 ;
01707       }
01708     
01709       // j suivants
01710       for (int j=nt-2 ; j > 0 ; j--) {
01711     // Positionnement
01712     xci -= 2*nr ;
01713     xco -= nr ;
01714     // Calcul
01715     for (int i=0 ; i<nr ; i++ ) {
01716       som[i] += 0.5*xci[i] ;
01717       xco[i] = som[i] ;
01718     }   // Fin de la boucle sur r
01719     xci += nr ;
01720     for (int i=0 ; i<nr ; i++ ) {
01721       som[i] = 0.5*xci[i] ;
01722     }
01723       }   // Fin de la boucle sur theta
01724       // j = 0 : sin(0*theta)
01725       xci -= nr ;
01726       xco -= nr ;
01727       for (int i=0; i<nr; i++) {
01728     xco[i] = 0. ;
01729       }
01730 
01731       // Positionnement phi suivant
01732       xci += nr*nt ;
01733       xco += nr*nt ;
01734       break;
01735     } // Fin du switch
01736   } // Fin de la boucle sur phi
01737   
01738   // On remet les choses la ou il faut
01739   delete [] tb->t ;
01740   tb->t = xo ;
01741   
01742   //Menage
01743   delete [] som ;
01744   
01745   // base de developpement
01746   int base_r = b & MSQ_R ;
01747   int base_p = b & MSQ_P ;
01748   switch(base_r){
01749   case(R_CHEBPI_P):
01750     b = R_CHEBPI_I | base_p | T_COSSIN_C ;
01751     break ;
01752   case(R_CHEBPI_I):
01753     b = R_CHEBPI_P | base_p | T_COSSIN_C ;
01754     break ;  
01755   default:
01756     b = base_r | base_p | T_COSSIN_C ;
01757     break;
01758   }
01759 }
01760 
01761             //--------------------- 
01762             // cas T_COSSIN_S
01763             //----------------
01764 void _mult_ct_t_cossin_s(Tbl* tb, int & b)
01765 {
01766   // Peut-etre rien a faire ?
01767   if (tb->get_etat() == ETATZERO) {
01768     int base_r = b & MSQ_R ;
01769     int base_p = b & MSQ_P ;
01770     switch(base_r){
01771     case(R_CHEBPI_P):
01772       b = R_CHEBPI_I | base_p | T_COSSIN_S ;
01773       break ;
01774     case(R_CHEBPI_I):
01775       b = R_CHEBPI_P | base_p | T_COSSIN_S ;
01776       break ;  
01777     default:
01778       b = base_r | base_p | T_COSSIN_S ;
01779       break;
01780     }
01781     return ;
01782   }
01783   
01784   // Protection
01785   assert(tb->get_etat() == ETATQCQ) ;
01786   
01787   // Pour le confort : nbre de points reels.
01788   int nr = (tb->dim).dim[0] ;
01789   int nt = (tb->dim).dim[1] ;
01790   int np = (tb->dim).dim[2] ;
01791   np = np - 2 ;
01792   
01793   // zone de travail
01794   double* som = new double [nr] ;
01795   
01796   // pt. sur le tableau de double resultat
01797   double* xo = new double[(tb->dim).taille] ;
01798   
01799     // Initialisation a zero :
01800     for (int i=0; i<(tb->dim).taille; i++) {
01801     xo[i] = 0 ; 
01802     }
01803     
01804   // On y va...
01805   double* xi = tb->t ;
01806   double* xci = xi ;    // Pointeurs
01807   double* xco = xo ;    //  courants
01808   
01809   // k = 0
01810   int m = 0 ;       // Parite de l'harmonique en phi
01811   // Dernier j: j = nt-1
01812   // Positionnement
01813   xci += nr * (nt-1) ;
01814   xco += nr * (nt-1) ;
01815     
01816   // Initialisation de som
01817   for (int i=0 ; i<nr ; i++) {
01818     som[i] = 0.5*xci[i] ;
01819     xco[i] = 0. ;
01820   }
01821     
01822   // j suivants
01823   for (int j=nt-2 ; j > 0 ; j--) {
01824     // Positionnement
01825     xci -= 2*nr ;
01826     xco -= nr ;
01827     // Calcul
01828     for (int i=0 ; i<nr ; i++ ) {
01829       som[i] += 0.5*xci[i] ;
01830       xco[i] = som[i] ;
01831     }   // Fin de la boucle sur r
01832     xci += nr ;
01833     for (int i=0 ; i<nr ; i++ ) {
01834        som[i] = 0.5*xci[i] ;
01835     }
01836   }   // Fin de la boucle sur theta
01837   // j = 0 : sin(0*theta)
01838   xci -= nr ;
01839   xco -= nr ;
01840   for (int i=0; i<nr; i++) {
01841     xco[i] = 0.0 ;
01842   }
01843   // Positionnement phi suivant
01844   xci += nr*nt ;
01845   xco += nr*nt ;
01846   
01847   // k=1 
01848   xci += nr*nt ;
01849   xco += nr*nt ;
01850   
01851   for (int k=2 ; k<np+1 ; k++) {
01852     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
01853     
01854     switch(m) {
01855     case 1:     // m impair: cos
01856       // Dernier j: j = nt-1
01857       // Positionnement
01858       xci += nr * (nt-1) ;
01859       xco += nr * (nt-1) ;
01860       
01861       // Initialisation de som
01862       for (int i=0 ; i<nr ; i++) {
01863     som[i] = 0.5*xci[i] ;
01864     xco[i] = 0.0 ;
01865       }
01866       
01867       // j suivants
01868       for (int j=nt-2 ; j > 0 ; j--) {
01869     // Positionnement
01870     xci -= 2*nr ;
01871     xco -= nr ;
01872     // Calcul
01873     for (int i=0 ; i<nr ; i++ ) {
01874       som[i] += 0.5*xci[i] ;
01875       xco[i] = som[i] ;
01876     }   // Fin de la boucle sur r
01877     xci += nr ;
01878     for (int i=0 ; i<nr ; i++ ) {
01879       som[i] = 0.5*xci[i] ;
01880     }
01881       }   // Fin de la boucle sur theta
01882       // premier theta
01883       // j=0 : le facteur 2...
01884       xci -= nr ;
01885       for (int i=0; i<nr; i++) {
01886       xco[i] += 0.5*xci[i] ;
01887       }
01888       xco -= nr ;
01889       for (int i=0 ; i<nr ; i++) {
01890     xco[i] = som[i] ;
01891       }
01892       // Positionnement phi suivant
01893       xci += nr*nt ;
01894       xco += nr*nt ;
01895       break ;
01896       
01897     case 0:     // m pair: sin
01898       // Dernier j: j = nt-1
01899       // Positionnement
01900       xci += nr * (nt-1) ;
01901       xco += nr * (nt-1) ;
01902       
01903       // Initialisation de som
01904       for (int i=0 ; i<nr ; i++) {
01905     som[i] = 0.5*xci[i] ;
01906     xco[i] = 0. ;
01907       }
01908       
01909       // j suivants
01910       for (int j=nt-2 ; j > 0 ; j--) {
01911     // Positionnement
01912     xci -= 2*nr ;
01913     xco -= nr ;
01914     // Calcul
01915     for (int i=0 ; i<nr ; i++ ) {
01916       som[i] += 0.5*xci[i] ;
01917       xco[i] = som[i] ;
01918     }   // Fin de la boucle sur r
01919     xci += nr ;
01920     for (int i=0 ; i<nr ; i++ ) {
01921       som[i] = 0.5*xci[i] ;
01922     }
01923       }   // Fin de la boucle sur theta
01924       // j = 0 : sin(0*theta)
01925       xci -= nr ;
01926       xco -= nr ;
01927       for (int i=0; i<nr; i++) {
01928     xco[i] = 0.0 ;
01929       }
01930       // Positionnement phi suivant
01931       xci += nr*nt ;
01932       xco += nr*nt ;
01933       break ;
01934     }   // Fin du switch
01935   }   // Fin de la boucle en phi
01936   
01937   // On remet les choses la ou il faut
01938   delete [] tb->t ;
01939   tb->t = xo ;
01940   
01941   //Menage
01942   delete [] som ;
01943   
01944   // base de developpement
01945   int base_r = b & MSQ_R ;
01946   int base_p = b & MSQ_P ;
01947   switch(base_r){
01948   case(R_CHEBPI_P):
01949     b = R_CHEBPI_I | base_p | T_COSSIN_S ;
01950     break ;
01951   case(R_CHEBPI_I):
01952     b = R_CHEBPI_P | base_p | T_COSSIN_S ;
01953     break ;  
01954   default:
01955     b = base_r | base_p | T_COSSIN_S ;
01956     break;
01957    }
01958 }
Generated on Tue Feb 7 01:35:18 2012 for LORENE by  doxygen 1.4.6