op_mult_st.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_st_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_st.C,v 1.6 2009/10/09 14:00:54 j_novak Exp $" ;
00024 
00025 /*
00026  * $Id: op_mult_st.C,v 1.6 2009/10/09 14:00:54 j_novak Exp $
00027  * $Log: op_mult_st.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:28  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:42:06  eric
00050  * Initialisation a zero du tableau xo avant le calcul.
00051  *
00052  * Revision 1.2  1999/11/23  16:13:53  novak
00053  * *** empty log message ***
00054  *
00055  * Revision 1.1  1999/11/23  14:31:50  novak
00056  * Initial revision
00057  *
00058  *
00059  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_st.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 sin theta
00065  * (Utilisation interne)
00066  * 
00067  *  void _mult_st_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_st_pas_prevu(Tbl * tb, int& base) {
00083   cout << "mult_st 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_st_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_SIN ;
00103       break ;
00104     case(R_CHEBPI_I):
00105       b = R_CHEBPI_P | base_p | T_SIN ;
00106       break ;  
00107     default:
00108       b = base_r | base_p | T_SIN ;
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     }
00166   }   // Fin de la boucle sur theta
00167       // j = 0 
00168   xci -= nr ;
00169   for (int i=0 ; i<nr ; i++ ) {
00170     xco[i] += 0.5*xci[i] ;
00171   }
00172   xco -= nr ;
00173   for (int i=0; i<nr; i++) {
00174     xco[i] = 0 ;
00175   }
00176   // Positionnement phi suivant
00177   xci += nr*nt ;
00178   xco += nr*nt ;
00179   
00180   // k = 1
00181   xci += nr*nt ;
00182   xco += nr*nt ;
00183   
00184   // k >= 2
00185   for (int k=2 ; k<np+1 ; k++) {
00186     
00187     // Dernier j: j = nt-1
00188     // Positionnement
00189     xci += nr * (nt-1) ;
00190     xco += nr * (nt-1) ;
00191     
00192     // Initialisation de som
00193     for (int i=0 ; i<nr ; i++) {
00194       som[i] = -0.5*xci[i] ;
00195       xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
00196     }
00197     
00198     // j suivants
00199     for (int j=nt-2 ; j > 0 ; j--) {
00200       // Positionnement
00201       xci -= 2*nr ;
00202       xco -= nr ;
00203       // Calcul
00204       for (int i=0 ; i<nr ; i++ ) {
00205     som[i] += 0.5*xci[i] ;
00206     xco[i] = som[i] ;
00207       } // Fin de la boucle sur r
00208       xci += nr ;
00209       for (int i=0; i<nr; i++) {
00210     som[i] = -0.5*xci[i] ;
00211       }
00212     }   // Fin de la boucle sur theta
00213     // j = 0
00214     xci -= nr ;
00215     for (int i = 0; i<nr; i++) {
00216       xco[i] += 0.5*xci[i] ;
00217     }
00218     xco -= nr ;
00219     for (int i=0; i<nr; i++) {
00220       xco[i] = 0 ;
00221     }
00222     // Positionnement phi suivant
00223     xci += nr*nt ;
00224     xco += nr*nt ;
00225   } // Fin de la boucle sur phi
00226   
00227   // On remet les choses la ou il faut
00228   delete [] tb->t ;
00229   tb->t = xo ;
00230   
00231   //Menage
00232   delete [] som ;
00233   
00234   // base de developpement
00235   int base_r = b & MSQ_R ;
00236   int base_p = b & MSQ_P ;
00237   switch(base_r){
00238   case(R_CHEBPI_P):
00239     b = R_CHEBPI_I | base_p | T_SIN ;
00240     break ;
00241   case(R_CHEBPI_I):
00242     b = R_CHEBPI_P | base_p | T_SIN ;
00243     break ;  
00244   default:
00245     b = base_r | base_p | T_SIN ;
00246       break;
00247   }
00248 }
00249             
00250             //------------
00251             // cas T_SIN
00252             //------------
00253             
00254 void _mult_st_t_sin(Tbl* tb, int & b)
00255 {
00256   // Peut-etre rien a faire ?
00257   if (tb->get_etat() == ETATZERO) {
00258     int base_r = b & MSQ_R ;
00259     int base_p = b & MSQ_P ;
00260     switch(base_r){
00261     case(R_CHEBPI_P):
00262       b = R_CHEBPI_I | base_p | T_COS ;
00263       break ;
00264     case(R_CHEBPI_I):
00265       b = R_CHEBPI_P | base_p | T_COS ;
00266       break ;  
00267     default:
00268       b = base_r | base_p | T_COS ;
00269       break;
00270     }
00271     return ;
00272   }
00273   
00274   // Protection
00275   assert(tb->get_etat() == ETATQCQ) ;
00276   
00277   // Pour le confort : nbre de points reels.
00278   int nr = (tb->dim).dim[0] ;
00279   int nt = (tb->dim).dim[1] ;
00280   int np = (tb->dim).dim[2] ;
00281   np = np - 2 ;
00282   
00283   // zone de travail
00284   double* som = new double [nr] ;
00285   
00286   // pt. sur le tableau de double resultat
00287   double* xo = new double[(tb->dim).taille] ;
00288   
00289     // Initialisation a zero :
00290     for (int i=0; i<(tb->dim).taille; i++) {
00291     xo[i] = 0 ; 
00292     }
00293     
00294   // On y va...
00295   double* xi = tb->t ;
00296   double* xci = xi ;    // Pointeurs
00297   double* xco = xo ;    //  courants
00298   
00299   // k = 0
00300   
00301   // Dernier j: j = nt-1
00302   // Positionnement
00303   xci += nr * (nt-1) ;
00304   xco += nr * (nt-1) ;
00305   
00306   // Initialisation de som
00307   for (int i=0 ; i<nr ; i++) {
00308     som[i] = 0.5*xci[i] ;
00309     xco[i] = 0. ;
00310   }
00311   
00312   // j suivants
00313   for (int j=nt-2 ; j > 0 ; j--) {
00314     // Positionnement
00315     xci -= 2*nr ;
00316     xco -= nr ;
00317     // Calcul
00318     for (int i=0 ; i<nr ; i++ ) {
00319       som[i] -= 0.5*xci[i] ;
00320       xco[i] = som[i] ;
00321     }   // Fin de la boucle sur r
00322     xci += nr ;
00323     for (int i=0; i<nr; i++) {
00324       som[i] = 0.5*xci[i] ;
00325     }
00326   }   // Fin de la boucle sur theta
00327       // j = 0
00328   xci -= nr ;
00329   xco -= nr  ;
00330   for (int i =0; i<nr ; i++) {
00331     xco[i] = som[i] ;
00332   }
00333   // Positionnement phi suivant
00334   xci += nr*nt ;
00335   xco += nr*nt ;
00336   
00337   // k = 1
00338   xci += nr*nt ;
00339   xco += nr*nt ;
00340   
00341   // k >= 2
00342   for (int k=2 ; k<np+1 ; k++) {
00343     
00344     // Dernier j: j = nt-1
00345     // Positionnement
00346     xci += nr * (nt-1) ;
00347     xco += nr * (nt-1) ;
00348     
00349     // Initialisation de som
00350     for (int i=0 ; i<nr ; i++) {
00351       som[i] = 0.5*xci[i] ;
00352       xco[i] = 0. ;
00353     }
00354     
00355     // j suivants
00356     for (int j=nt-2 ; j > 0 ; j--) {
00357       // Positionnement
00358       xci -= 2*nr ;
00359       xco -= nr ;
00360       // Calcul
00361       for (int i=0 ; i<nr ; i++ ) {
00362     som[i] -= 0.5*xci[i] ;
00363     xco[i] = som[i] ;
00364       } // Fin de la boucle sur r
00365     xci += nr ;
00366     for (int i=0 ; i<nr ; i++ ) {
00367       som[i] = 0.5*xci[i] ;
00368     }
00369     }   // Fin de la boucle sur theta
00370     // j = 0 
00371     xci -= nr ;
00372     xco -= nr ;
00373     for (int i=0; i<nr; i++) {
00374       xco[i] = som[i] ;
00375     }
00376     // Positionnement phi suivant
00377     xci += nr*nt ;
00378     xco += nr*nt ;
00379   } // Fin de la boucle sur phi
00380   
00381   // On remet les choses la ou il faut
00382   delete [] tb->t ;
00383   tb->t = xo ;
00384   
00385   //Menage
00386   delete [] som ;
00387   
00388   // base de developpement
00389   int base_r = b & MSQ_R ;
00390   int base_p = b & MSQ_P ;
00391   switch(base_r){
00392   case(R_CHEBPI_P):
00393     b = R_CHEBPI_I | base_p | T_COS ;
00394     break ;
00395   case(R_CHEBPI_I):
00396     b = R_CHEBPI_P | base_p | T_COS ;
00397     break ;  
00398   default:
00399     b = base_r | base_p | T_COS ;
00400     break;
00401   }
00402 }
00403             //----------------
00404             // cas T_COS_P
00405             //----------------
00406             
00407 void _mult_st_t_cos_p(Tbl* tb, int & b)
00408 {
00409 
00410   // Peut-etre rien a faire ?
00411   if (tb->get_etat() == ETATZERO) {
00412     int base_r = b & MSQ_R ;
00413     int base_p = b & MSQ_P ;
00414     b = base_r | base_p | T_SIN_I ;
00415     return ;
00416   }
00417     
00418   // Protection
00419   assert(tb->get_etat() == ETATQCQ) ;
00420   
00421   // Pour le confort : nbre de points reels.
00422   int nr = (tb->dim).dim[0] ;
00423   int nt = (tb->dim).dim[1] ;
00424   int np = (tb->dim).dim[2] ;
00425   np = np - 2 ;
00426   
00427   // zone de travail
00428   double* som = new double [nr] ;
00429   
00430   // pt. sur le tableau de double resultat
00431   double* xo = new double[(tb->dim).taille] ;
00432   
00433     // Initialisation a zero :
00434     for (int i=0; i<(tb->dim).taille; i++) {
00435     xo[i] = 0 ; 
00436     }
00437     
00438   // On y va...
00439   double* xi = tb->t ;
00440   double* xci = xi ;    // Pointeurs
00441   double* xco = xo ;    //  courants
00442   
00443   // k = 0
00444   
00445   // Dernier j: j = nt-1
00446   // Positionnement
00447   xci += nr * (nt-1) ;
00448   xco += nr * (nt-1) ;
00449   
00450   // Initialisation de som 
00451   for (int i=0 ; i<nr ; i++) {
00452     som[i] = -0.5*xci[i] ;
00453     xco[i] = 0. ;   // mise a zero du dernier coef en theta
00454   }
00455   
00456   // j suivants
00457   for (int j=nt-2 ; j > 0 ; j--) {
00458     // Positionnement
00459     xci -= nr ;
00460     xco -= nr ;
00461     // Calcul
00462     for (int i=0 ; i<nr ; i++ ) {
00463       som[i] += 0.5*xci[i] ;
00464       xco[i] = som[i] ;
00465       som[i] = -0.5*xci[i] ;
00466     }   // Fin de la boucle sur r
00467   }   // Fin de la boucle sur theta
00468   if (nt > 1) {
00469       // j = 0 
00470       xci -= nr ;
00471       xco -= nr ;
00472       for (int i=0 ; i<nr ; i++ ) {
00473       xco[i] = som[i]+xci[i] ;
00474       } // Fin de la boucle sur r
00475   }
00476   // Positionnement phi suivant
00477   xci += nr*nt ;
00478   xco += nr*nt ;
00479   
00480   // k = 1
00481   xci += nr*nt ;
00482   xco += nr*nt ;
00483   
00484   // k >= 2
00485   for (int k=2 ; k<np+1 ; k++) {
00486     
00487     // Dernier j: j = nt-1
00488     // Positionnement
00489     xci += nr * (nt-1) ;
00490     xco += nr * (nt-1) ;
00491     
00492     // Initialisation de som
00493     for (int i=0 ; i<nr ; i++) {
00494       som[i] = -0.5*xci[i] ;
00495       xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
00496     }
00497     
00498     // j suivants
00499     for (int j=nt-2 ; j > 0 ; j--) {
00500       // Positionnement
00501       xci -= nr ;
00502       xco -= nr ;
00503       // Calcul
00504       for (int i=0 ; i<nr ; i++ ) {
00505     som[i] += 0.5*xci[i] ;
00506     xco[i] = som[i] ;
00507     som[i] = -0.5*xci[i] ;
00508       } // Fin de la boucle sur r
00509     }   // Fin de la boucle sur theta
00510     // j = 0
00511     xci -= nr ;
00512     xco -= nr ;
00513     for (int i = 0; i<nr; i++) {
00514       xco[i] = xci[i] + som[i] ;
00515     }
00516     // Positionnement phi suivant
00517     xci += nr*nt ;
00518     xco += nr*nt ;
00519   } // Fin de la boucle sur phi
00520   
00521   // On remet les choses la ou il faut
00522   delete [] tb->t ;
00523   tb->t = xo ;
00524   
00525   //Menage
00526   delete [] som ;
00527   
00528   // base de developpement
00529   int base_r = b & MSQ_R ;
00530   int base_p = b & MSQ_P ;
00531   b = base_r | base_p | T_SIN_I ;
00532 }
00533             
00534             //--------------
00535             // cas T_SIN_P
00536             //--------------
00537             
00538 void _mult_st_t_sin_p(Tbl* tb, int & b)
00539 {
00540   // Peut-etre rien a faire ?
00541   if (tb->get_etat() == ETATZERO) {
00542     int base_r = b & MSQ_R ;
00543     int base_p = b & MSQ_P ;
00544     b = base_r | base_p | T_COS_I ;
00545     return ;
00546   }
00547   
00548   // Protection
00549   assert(tb->get_etat() == ETATQCQ) ;
00550   
00551   // Pour le confort : nbre de points reels.
00552   int nr = (tb->dim).dim[0] ;
00553   int nt = (tb->dim).dim[1] ;
00554   int np = (tb->dim).dim[2] ;
00555   np = np - 2 ;
00556   
00557   // zone de travail
00558   double* som = new double [nr] ;
00559   
00560   // pt. sur le tableau de double resultat
00561   double* xo = new double[(tb->dim).taille] ;
00562   
00563     // Initialisation a zero :
00564     for (int i=0; i<(tb->dim).taille; i++) {
00565     xo[i] = 0 ; 
00566     }
00567     
00568   // On y va...
00569   double* xi = tb->t ;
00570   double* xci = xi ;    // Pointeurs
00571   double* xco = xo ;    //  courants
00572   
00573   // k = 0
00574   
00575   // Dernier j: j = nt-1
00576   // Positionnement
00577   xci += nr * (nt-1) ;
00578   xco += nr * (nt-1) ;
00579   
00580   // Initialisation de som
00581   for (int i=0 ; i<nr ; i++) {
00582     som[i] = 0.5*xci[i] ;
00583     xco[i] = 0. ;
00584   }
00585   
00586   // j suivants
00587   for (int j=nt-2 ; j > 0 ; j--) {
00588     // Positionnement
00589     xci -= nr ;
00590     xco -= nr ;
00591     // Calcul
00592     for (int i=0 ; i<nr ; i++ ) {
00593       som[i] -= 0.5*xci[i] ;
00594       xco[i] = som[i] ;
00595       som[i] = 0.5*xci[i] ;
00596     }   // Fin de la boucle sur r
00597   }   // Fin de la boucle sur theta
00598   if (nt > 1) {
00599       // j = 0
00600       xci -= nr ;
00601       xco -= nr  ;
00602       for (int i =0; i<nr ; i++) {
00603       xco[i] = som[i] ;
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_COS_I ;
00662   
00663 }
00664             //--------------
00665             // cas T_SIN_I
00666             //--------------
00667             
00668 void _mult_st_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_COS_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   // premier theta
00728   for (int i=0 ; i<nr ; i++) {
00729     xco[i] = som[i] ;
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     // premier theta
00765     for (int i=0 ; i<nr ; i++) {
00766       xco[i] = som[i] ;
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_COS_P ;
00784   
00785 }
00786                     //---------------------
00787                     // cas T_COS_I
00788             //---------------------
00789 void _mult_st_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_SIN_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   for (int i=0; i<nr; i++) {
00849     xco[i] = 0. ;
00850   }
00851   // Positionnement phi suivant
00852   xci += nr*nt ;
00853   xco += nr*nt ;
00854   
00855   // k = 1
00856   xci += nr*nt ;
00857   xco += nr*nt ;
00858     
00859   // k >= 2
00860   for (int k=2 ; k<np+1 ; k++) {
00861     
00862     // Dernier j: j = nt-1
00863     // Positionnement
00864     xci += nr * (nt-1) ;
00865     xco += nr * (nt-1) ;
00866     
00867     // Initialisation de som
00868     for (int i=0 ; i<nr ; i++) {
00869       som[i] = 0. ;
00870     }
00871     
00872     // j suivants
00873     for (int j=nt-1 ; j > 0 ; j--) {
00874       // Positionnement
00875       xci -= nr ;
00876       // Calcul
00877       for (int i=0 ; i<nr ; i++ ) {
00878     som[i] += 0.5*xci[i] ;
00879     xco[i] = som[i] ;
00880     som[i] = -0.5*xci[i] ;
00881       } // Fin de la boucle sur r
00882       xco -= nr ;
00883     }   // Fin de la boucle sur theta
00884   for (int i=0; i<nr; i++) {
00885     xco[i] = 0. ;
00886   }
00887 
00888     // Positionnement phi suivant
00889     xci += nr*nt ;
00890     xco += nr*nt ;
00891   } // Fin de la boucle sur phi
00892   
00893   // On remet les choses la ou il faut
00894   delete [] tb->t ;
00895   tb->t = xo ;
00896   
00897   //Menage
00898   delete [] som ;
00899   
00900   // base de developpement
00901   int base_r = b & MSQ_R ;
00902   int base_p = b & MSQ_P ;
00903   b = base_r | base_p | T_SIN_P ;
00904   
00905 }
00906             //----------------------
00907                         // cas T_COSSIN_CP
00908             //----------------------
00909 void _mult_st_t_cossin_cp(Tbl* tb, int & b)
00910 {
00911   // Peut-etre rien a faire ?
00912   if (tb->get_etat() == ETATZERO) {
00913     int base_r = b & MSQ_R ;
00914     int base_p = b & MSQ_P ;
00915     b = base_r | base_p | T_COSSIN_SI ;
00916     return ;
00917   }
00918   
00919   // Protection
00920   assert(tb->get_etat() == ETATQCQ) ;
00921   
00922   // Pour le confort : nbre de points reels.
00923   int nr = (tb->dim).dim[0] ;
00924   int nt = (tb->dim).dim[1] ;
00925   int np = (tb->dim).dim[2] ;
00926   np = np - 2 ;
00927     
00928   // zone de travail
00929   double* som = new double [nr] ;
00930   
00931   // pt. sur le tableau de double resultat
00932   double* xo = new double[(tb->dim).taille] ;
00933     
00934     // Initialisation a zero :
00935     for (int i=0; i<(tb->dim).taille; i++) {
00936     xo[i] = 0 ; 
00937     }
00938     
00939   // On y va...
00940   double* xi = tb->t ;
00941   double* xci = xi ;    // Pointeurs
00942   double* xco = xo ;    //  courants
00943   
00944   // k = 0
00945   int m = 0 ;       // Parite de l'harmonique en phi
00946   // Dernier j: j = nt-1
00947   // Positionnement
00948   xci += nr * (nt-1) ;
00949   xco += nr * (nt-1) ;
00950     
00951   // Initialisation de som
00952   for (int i=0 ; i<nr ; i++) {
00953     som[i] = -0.5*xci[i] ;
00954     xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
00955   }
00956     
00957   // j suivants
00958   for (int j=nt-2 ; j > 0 ; j--) {
00959     // Positionnement
00960     xci -= nr ;
00961     xco -= nr ;
00962     // Calcul
00963     for (int i=0 ; i<nr ; i++ ) {
00964       som[i] += 0.5*xci[i] ;
00965       xco[i] = som[i] ;
00966       som[i] = -0.5*xci[i] ;
00967     }   // Fin de la boucle sur r
00968   }   // Fin de la boucle sur theta
00969   // j = 0
00970   xci -= nr ;
00971   xco -= nr ;
00972   for (int i = 0; i<nr; i++) {
00973     xco[i] = xci[i] + som[i] ;
00974   }
00975   // Positionnement phi suivant
00976   xci += nr*nt ;
00977   xco += nr*nt ;
00978   
00979   // k=1
00980   xci += nr*nt ;
00981   xco += nr*nt ;
00982   
00983   // k >= 2
00984   for (int k=2 ; k<np+1 ; k++) {
00985     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
00986     
00987     switch(m) {
00988     case 0:     // m pair: cos(pair)
00989       // Dernier j: j = nt-1
00990       // Positionnement
00991       xci += nr * (nt-1) ;
00992       xco += nr * (nt-1) ;
00993       
00994       // Initialisation de som
00995       for (int i=0 ; i<nr ; i++) {
00996     som[i] = -0.5*xci[i] ;
00997     xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
00998       }
00999       
01000       // j suivants
01001       for (int j=nt-2 ; j > 0 ; j--) {
01002     // Positionnement
01003     xci -= nr ;
01004     xco -= nr ;
01005     // Calcul
01006     for (int i=0 ; i<nr ; i++ ) {
01007       som[i] += 0.5*xci[i] ;
01008       xco[i] = som[i] ;
01009       som[i] = -0.5*xci[i] ;
01010     }   // Fin de la boucle sur r
01011       }   // Fin de la boucle sur theta
01012       // j = 0
01013       xci -= nr ;
01014       xco -= nr ;
01015       for (int i = 0; i<nr; i++) {
01016     xco[i] = xci[i] + som[i] ;
01017       }
01018       // Positionnement phi suivant
01019       xci += nr*nt ;
01020       xco += nr*nt ;
01021       break ;
01022       
01023     case 1:     // m impair: sin(impair)
01024       // Dernier j: j = nt-1
01025       // Positionnement
01026       xci += nr * (nt-1) ;
01027       xco += nr * (nt-1) ;
01028       
01029       // Initialisation de som
01030       for (int i=0 ; i<nr ; i++) {
01031     som[i] = 0. ;
01032       }
01033     
01034       // j suivants
01035       for (int j=nt-1 ; j > 0 ; j--) {
01036     // Positionnement
01037     xci -= nr ;
01038     // Calcul
01039     for (int i=0 ; i<nr ; i++ ) {
01040       som[i] -= 0.5*xci[i] ;
01041       xco[i] = som[i] ;
01042       som[i] = 0.5*xci[i] ;
01043     }   // Fin de la boucle sur r
01044     xco -= nr ;
01045       }   // Fin de la boucle sur theta
01046       // premier theta
01047       for (int i=0 ; i<nr ; i++) {
01048     xco[i] = som[i] ;
01049       }
01050       // Positionnement phi suivant
01051       xci += nr*nt ;
01052       xco += nr*nt ;
01053       break;
01054     } // Fin du switch
01055   } // Fin de la boucle sur phi
01056   
01057   // On remet les choses la ou il faut
01058   delete [] tb->t ;
01059   tb->t = xo ;
01060   
01061   //Menage
01062   delete [] som ;
01063   
01064   // base de developpement
01065   int base_r = b & MSQ_R ;
01066   int base_p = b & MSQ_P ;
01067   b = base_r | base_p | T_COSSIN_SI ;
01068 }
01069             
01070             //------------------
01071             // cas T_COSSIN_CI
01072             //----------------
01073 void _mult_st_t_cossin_ci(Tbl* tb, int & b)
01074 {
01075   // Peut-etre rien a faire ?
01076   if (tb->get_etat() == ETATZERO) {
01077     int base_r = b & MSQ_R ;
01078     int base_p = b & MSQ_P ;
01079     b = base_r | base_p | T_COSSIN_SP ;
01080     return ;
01081   }
01082   
01083   // Protection
01084   assert(tb->get_etat() == ETATQCQ) ;
01085   
01086   // Pour le confort : nbre de points reels.
01087   int nr = (tb->dim).dim[0] ;
01088   int nt = (tb->dim).dim[1] ;
01089   int np = (tb->dim).dim[2] ;
01090   np = np - 2 ;
01091   
01092   // zone de travail
01093   double* som = new double [nr] ;
01094   
01095   // pt. sur le tableau de double resultat
01096   double* xo = new double[(tb->dim).taille] ;
01097   
01098     // Initialisation a zero :
01099     for (int i=0; i<(tb->dim).taille; i++) {
01100     xo[i] = 0 ; 
01101     }
01102     
01103   // On y va...
01104   double* xi = tb->t ;
01105   double* xci = xi ;    // Pointeurs
01106   double* xco = xo ;    //  courants
01107   
01108   // k = 0
01109   int m = 0 ;       // Parite de l'harmonique en phi
01110   // Dernier j: j = nt-1
01111   // Positionnement
01112   xci += nr * (nt-1) ;
01113   xco += nr * (nt-1) ;
01114   
01115   // Initialisation de som
01116   for (int i=0 ; i<nr ; i++) {
01117     som[i] = 0. ;
01118   }
01119   
01120   // j suivants
01121   for (int j=nt-1 ; j > 0 ; j--) {
01122     // Positionnement
01123     xci -= nr ;
01124     // Calcul
01125     for (int i=0 ; i<nr ; i++ ) {
01126       som[i] += 0.5*xci[i] ;
01127       xco[i] = som[i] ;
01128       som[i] = -0.5*xci[i] ;
01129     }   // Fin de la boucle sur r
01130     xco -= nr ;
01131   }   // Fin de la boucle sur theta
01132   for (int i=0; i<nr; i++) {
01133     xco[i] = 0. ;
01134   }
01135 
01136   // Positionnement phi suivant
01137   xci += nr*nt ;
01138   xco += nr*nt ;
01139   
01140   // k=1
01141   xci += nr*nt ;
01142   xco += nr*nt ;
01143   
01144   // k >= 2
01145   for (int k=2 ; k<np+1 ; k++) {
01146     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
01147     
01148     switch(m) {
01149     case 0:     // m pair: cos(impair)
01150       // Dernier j: j = nt-1
01151       // Positionnement
01152       xci += nr * (nt-1) ;
01153       xco += nr * (nt-1) ;
01154       
01155       // Initialisation de som
01156       for (int i=0 ; i<nr ; i++) {
01157     som[i] = 0. ;
01158       }
01159       
01160       // j suivants
01161       for (int j=nt-1 ; j > 0 ; j--) {
01162     // Positionnement
01163     xci -= nr ;
01164     // Calcul
01165     for (int i=0 ; i<nr ; i++ ) {
01166       som[i] += 0.5*xci[i] ;
01167       xco[i] = som[i] ;
01168       som[i] = -0.5*xci[i] ;
01169     }   // Fin de la boucle sur r
01170     xco -= nr ;
01171       }   // Fin de la boucle sur theta
01172   for (int i=0; i<nr; i++) {
01173     xco[i] = 0. ;
01174   }
01175 
01176       // Positionnement phi suivant
01177       xci += nr*nt ;
01178       xco += nr*nt ;
01179       break ;
01180       
01181     case 1:
01182       // Dernier j: j = nt-1
01183       // Positionnement
01184       xci += nr * (nt-1) ;
01185       xco += nr * (nt-1) ;
01186       
01187       // Initialisation de som
01188       for (int i=0 ; i<nr ; i++) {
01189     som[i] = 0.5*xci[i] ;
01190     xco[i] = 0. ;
01191       }
01192       
01193       // j suivants
01194       for (int j=nt-2 ; j > 0 ; j--) {
01195     // Positionnement
01196     xci -= nr ;
01197     xco -= nr ;
01198     // Calcul
01199     for (int i=0 ; i<nr ; i++ ) {
01200       som[i] -= 0.5*xci[i] ;
01201       xco[i] = som[i] ;
01202       som[i] = 0.5*xci[i] ;
01203     }   // Fin de la boucle sur r
01204       }   // Fin de la boucle sur theta
01205       // j = 0 
01206       xci -= nr ;
01207       xco -= nr ;
01208       for (int i=0; i<nr; i++) {
01209     xco[i] = som[i] ;
01210       }
01211       // Positionnement phi suivant
01212       xci += nr*nt ;
01213       xco += nr*nt ;
01214       break ;
01215     }   // Fin du switch
01216   }   // Fin de la boucle en phi
01217   
01218   // On remet les choses la ou il faut
01219   delete [] tb->t ;
01220   tb->t = xo ;
01221   
01222   //Menage
01223   delete [] som ;
01224   
01225   // base de developpement
01226   int base_r = b & MSQ_R ;
01227   int base_p = b & MSQ_P ;
01228   b = base_r | base_p | T_COSSIN_SP ;
01229 }
01230 
01231             //--------------------- 
01232             // cas T_COSSIN_SI
01233             //----------------
01234 void _mult_st_t_cossin_si(Tbl* tb, int & b)
01235 {
01236   // Peut-etre rien a faire ?
01237   if (tb->get_etat() == ETATZERO) {
01238     int base_r = b & MSQ_R ;
01239     int base_p = b & MSQ_P ;
01240     b = base_r | base_p | T_COSSIN_CP ;
01241     return ;
01242   }
01243   
01244   // Protection
01245   assert(tb->get_etat() == ETATQCQ) ;
01246   
01247   // Pour le confort : nbre de points reels.
01248   int nr = (tb->dim).dim[0] ;
01249   int nt = (tb->dim).dim[1] ;
01250   int np = (tb->dim).dim[2] ;
01251   np = np - 2 ;
01252   
01253   // zone de travail
01254   double* som = new double [nr] ;
01255   
01256   // pt. sur le tableau de double resultat
01257   double* xo = new double[(tb->dim).taille] ;
01258   
01259     // Initialisation a zero :
01260     for (int i=0; i<(tb->dim).taille; i++) {
01261     xo[i] = 0 ; 
01262     }
01263     
01264   // On y va...
01265   double* xi = tb->t ;
01266   double* xci = xi ;    // Pointeurs
01267   double* xco = xo ;    //  courants
01268   
01269   // k = 0
01270   int m = 0 ;       // Parite de l'harmonique en phi
01271   // Dernier j: j = nt-1
01272   // Positionnement
01273   xci += nr * (nt-1) ;
01274   xco += nr * (nt-1) ;
01275   
01276   // Initialisation de som
01277   for (int i=0 ; i<nr ; i++) {
01278     som[i] = 0. ;
01279   }
01280   
01281   // j suivants
01282   for (int j=nt-1 ; j > 0 ; j--) {
01283     // Positionnement
01284     xci -= nr ;
01285     // Calcul
01286     for (int i=0 ; i<nr ; i++ ) {
01287       som[i] -= 0.5*xci[i] ;
01288       xco[i] = som[i] ;
01289       som[i] = 0.5*xci[i] ;
01290     }   // Fin de la boucle sur r
01291     xco -= nr ;
01292   }   // Fin de la boucle sur theta
01293   // premier theta
01294   for (int i=0 ; i<nr ; i++) {
01295     xco[i] = som[i] ;
01296   }
01297   // Positionnement phi suivant
01298   xci += nr*nt ;
01299   xco += nr*nt ;
01300   
01301   // k=1 
01302   xci += nr*nt ;
01303   xco += nr*nt ;
01304   
01305   // k >= 2
01306   for (int k=2 ; k<np+1 ; k++) {
01307     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
01308     
01309     switch(m) {
01310     case 0:     // m pair: sin(impair)
01311       // Dernier j: j = nt-1
01312       // Positionnement
01313       xci += nr * (nt-1) ;
01314       xco += nr * (nt-1) ;
01315     
01316       // Initialisation de som
01317       for (int i=0 ; i<nr ; i++) {
01318     som[i] = 0. ;
01319       }
01320       
01321       // j suivants
01322       for (int j=nt-1 ; j > 0 ; j--) {
01323     // Positionnement
01324     xci -= nr ;
01325         // Calcul
01326     for (int i=0 ; i<nr ; i++ ) {
01327       som[i] -= 0.5*xci[i] ;
01328       xco[i] = som[i] ;
01329       som[i] = 0.5*xci[i] ;
01330     }   // Fin de la boucle sur r
01331     xco -= nr ;
01332       }   // Fin de la boucle sur theta
01333       // premier theta
01334       for (int i=0 ; i<nr ; i++) {
01335     xco[i] = som[i] ;
01336       }
01337       // Positionnement phi suivant
01338       xci += nr*nt ;
01339       xco += nr*nt ;
01340       break ;
01341       
01342     case 1: // m impair cos(pair)
01343       // Dernier j: j = nt-1
01344       // Positionnement
01345       xci += nr * (nt-1) ;
01346       xco += nr * (nt-1) ;
01347       
01348       // Initialisation de som
01349       for (int i=0 ; i<nr ; i++) {
01350     som[i] = -0.5*xci[i] ;
01351     xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
01352       }
01353       
01354       // j suivants
01355       for (int j=nt-2 ; j > 0 ; j--) {
01356     // Positionnement
01357     xci -= nr ;
01358     xco -= nr ;
01359     // Calcul
01360     for (int i=0 ; i<nr ; i++ ) {
01361       som[i] += 0.5*xci[i] ;
01362       xco[i] = som[i] ;
01363       som[i] = -0.5*xci[i] ;
01364     }   // Fin de la boucle sur r
01365       }   // Fin de la boucle sur theta
01366       // j = 0
01367       xci -= nr ;
01368       xco -= nr ;
01369       for (int i = 0; i<nr; i++) {
01370     xco[i] = xci[i] + som[i] ;
01371       }
01372       // Positionnement phi suivant
01373       xci += nr*nt ;
01374       xco += nr*nt ;
01375       break ;
01376     }   // Fin du switch
01377   }   // Fin de la boucle en phi
01378   
01379   // On remet les choses la ou il faut
01380   delete [] tb->t ;
01381   tb->t = xo ;
01382   
01383   //Menage
01384   delete [] som ;
01385   
01386   // base de developpement
01387   int base_r = b & MSQ_R ;
01388   int base_p = b & MSQ_P ;
01389   b = base_r | base_p | T_COSSIN_CP ;
01390   
01391   
01392 }
01393             //--------------------- 
01394             // cas T_COSSIN_SP
01395             //----------------
01396 void _mult_st_t_cossin_sp(Tbl* tb, int & b)
01397 {
01398   // Peut-etre rien a faire ?
01399   if (tb->get_etat() == ETATZERO) {
01400     int base_r = b & MSQ_R ;
01401     int base_p = b & MSQ_P ;
01402     b = base_r | base_p | T_COSSIN_CI ;
01403     return ;
01404   }
01405   
01406   // Protection
01407   assert(tb->get_etat() == ETATQCQ) ;
01408   
01409   // Pour le confort : nbre de points reels.
01410   int nr = (tb->dim).dim[0] ;
01411   int nt = (tb->dim).dim[1] ;
01412   int np = (tb->dim).dim[2] ;
01413   np = np - 2 ;
01414   
01415   // zone de travail
01416   double* som = new double [nr] ;
01417   
01418   // pt. sur le tableau de double resultat
01419   double* xo = new double[(tb->dim).taille] ;
01420   
01421     // Initialisation a zero :
01422     for (int i=0; i<(tb->dim).taille; i++) {
01423     xo[i] = 0 ; 
01424     }
01425     
01426   // On y va...
01427   double* xi = tb->t ;
01428   double* xci = xi ;    // Pointeurs
01429   double* xco = xo ;    //  courants
01430   
01431   // k = 0
01432   int m = 0 ;       // Parite de l'harmonique en phi
01433   // Dernier j: j = nt-1
01434   // Positionnement
01435   xci += nr * (nt-1) ;
01436   xco += nr * (nt-1) ;
01437     
01438   // Initialisation de som
01439   for (int i=0 ; i<nr ; i++) {
01440     som[i] = 0.5*xci[i] ;
01441     xco[i] = 0. ;
01442   }
01443     
01444   // j suivants
01445   for (int j=nt-2 ; j > 0 ; j--) {
01446     // Positionnement
01447     xci -= nr ;
01448     xco -= nr ;
01449     // Calcul
01450     for (int i=0 ; i<nr ; i++ ) {
01451       som[i] -= 0.5*xci[i] ;
01452       xco[i] = som[i] ;
01453       som[i] = 0.5*xci[i] ;
01454     }   // Fin de la boucle sur r
01455   }   // Fin de la boucle sur theta
01456   // j = 0 
01457   xci -= nr ;
01458   xco -= nr ;
01459   for (int i=0; i<nr; i++) {
01460     xco[i] = som[i] ;
01461   }
01462   // Positionnement phi suivant
01463   xci += nr*nt ;
01464   xco += nr*nt ;
01465   
01466   // k=1 
01467   xci += nr*nt ;
01468   xco += nr*nt ;
01469   
01470   for (int k=2 ; k<np+1 ; k++) {
01471     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
01472     
01473     switch(m) {
01474     case 1:     // m impair: cos(impair)
01475       // Dernier j: j = nt-1
01476       // Positionnement
01477       xci += nr * (nt-1) ;
01478       xco += nr * (nt-1) ;
01479       
01480       // Initialisation de som
01481       for (int i=0 ; i<nr ; i++) {
01482     som[i] = 0. ;
01483       }
01484       
01485       // j suivants
01486       for (int j=nt-1 ; j > 0 ; j--) {
01487     // Positionnement
01488     xci -= nr ;
01489     // Calcul
01490     for (int i=0 ; i<nr ; i++ ) {
01491       som[i] += 0.5*xci[i] ;
01492       xco[i] = som[i] ;
01493       som[i] = -0.5*xci[i] ;
01494     }   // Fin de la boucle sur r
01495     xco -= nr ;
01496       }   // Fin de la boucle sur theta
01497   for (int i=0; i<nr; i++) {
01498     xco[i] = 0. ;
01499   }
01500 
01501       // Positionnement phi suivant
01502       xci += nr*nt ;
01503       xco += nr*nt ;
01504       break ;
01505       
01506     case 0:     // m pair: sin(pair)
01507       // Dernier j: j = nt-1
01508       // Positionnement
01509       xci += nr * (nt-1) ;
01510       xco += nr * (nt-1) ;
01511       
01512       // Initialisation de som
01513       for (int i=0 ; i<nr ; i++) {
01514     som[i] = 0.5*xci[i] ;
01515     xco[i] = 0. ;
01516       }
01517       
01518       // j suivants
01519       for (int j=nt-2 ; j > 0 ; j--) {
01520     // Positionnement
01521     xci -= nr ;
01522     xco -= nr ;
01523     // Calcul
01524     for (int i=0 ; i<nr ; i++ ) {
01525       som[i] -= 0.5*xci[i] ;
01526       xco[i] = som[i] ;
01527       som[i] = 0.5*xci[i] ;
01528     }   // Fin de la boucle sur r
01529       }   // Fin de la boucle sur theta
01530       // j = 0 
01531       xci -= nr ;
01532       xco -= nr ;
01533       for (int i=0; i<nr; i++) {
01534     xco[i] = som[i] ;
01535       }
01536       // Positionnement phi suivant
01537       xci += nr*nt ;
01538       xco += nr*nt ;
01539       break ;
01540     }   // Fin du switch
01541   }   // Fin de la boucle en phi
01542   
01543   // On remet les choses la ou il faut
01544   delete [] tb->t ;
01545   tb->t = xo ;
01546   
01547   //Menage
01548   delete [] som ;
01549   
01550   // base de developpement
01551   int base_r = b & MSQ_R ;
01552   int base_p = b & MSQ_P ;
01553   b = base_r | base_p | T_COSSIN_CI ;
01554     
01555 }
01556 
01557             //----------------------
01558                         // cas T_COSSIN_C
01559             //----------------------
01560 void _mult_st_t_cossin_c(Tbl* tb, int & b)
01561 {
01562   // Peut-etre rien a faire ?
01563   if (tb->get_etat() == ETATZERO) {
01564     int base_r = b & MSQ_R ;
01565     int base_p = b & MSQ_P ;
01566     switch(base_r){
01567     case(R_CHEBPI_P):
01568       b = R_CHEBPI_I | base_p | T_COSSIN_S ;
01569       break ;
01570     case(R_CHEBPI_I):
01571       b = R_CHEBPI_P | base_p | T_COSSIN_S ;
01572       break ;  
01573     default:
01574       b = base_r | base_p | T_COSSIN_S ;
01575       break;
01576     }
01577     return ;
01578   }
01579   
01580   // Protection
01581   assert(tb->get_etat() == ETATQCQ) ;
01582   
01583   // Pour le confort : nbre de points reels.
01584   int nr = (tb->dim).dim[0] ;
01585   int nt = (tb->dim).dim[1] ;
01586   int np = (tb->dim).dim[2] ;
01587   np = np - 2 ;
01588     
01589   // zone de travail
01590   double* som = new double [nr] ;
01591   
01592   // pt. sur le tableau de double resultat
01593   double* xo = new double[(tb->dim).taille] ;
01594     
01595     // Initialisation a zero :
01596     for (int i=0; i<(tb->dim).taille; i++) {
01597     xo[i] = 0 ; 
01598     }
01599     
01600   // On y va...
01601   double* xi = tb->t ;
01602   double* xci = xi ;    // Pointeurs
01603   double* xco = xo ;    //  courants
01604   
01605   // k = 0
01606   int m = 0 ;       // Parite de l'harmonique en phi
01607   // Dernier j: j = nt-1
01608   // Positionnement
01609   xci += nr * (nt-1) ;
01610   xco += nr * (nt-1) ;
01611     
01612   // Initialisation de som
01613   for (int i=0 ; i<nr ; i++) {
01614     som[i] = -0.5*xci[i] ;
01615     xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
01616   }
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
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] = 0. ;
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
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] = 0. ;
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. ; // mise a zero dui dernier coefficient en theta.
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       xci -= nr;
01725       xco -= nr;
01726       // premier theta
01727       for (int i=0 ; i<nr ; i++) {
01728     xco[i] = som[i] ;
01729       }
01730       // Positionnement phi suivant
01731       xci += nr*nt ;
01732       xco += nr*nt ;
01733       break;
01734     } // Fin du switch
01735   } // Fin de la boucle sur phi
01736   
01737   // On remet les choses la ou il faut
01738   delete [] tb->t ;
01739   tb->t = xo ;
01740   
01741   //Menage
01742   delete [] som ;
01743   
01744   // base de developpement
01745   int base_r = b & MSQ_R ;
01746   int base_p = b & MSQ_P ;
01747   switch(base_r){
01748   case(R_CHEBPI_P):
01749     b = R_CHEBPI_I | base_p | T_COSSIN_S ;
01750     break ;
01751   case(R_CHEBPI_I):
01752     b = R_CHEBPI_P | base_p | T_COSSIN_S ;
01753     break ;  
01754   default:
01755     b = base_r | base_p | T_COSSIN_S ;
01756     break;
01757   }
01758 }
01759 
01760             //--------------------- 
01761             // cas T_COSSIN_S
01762             //----------------
01763 void _mult_st_t_cossin_s(Tbl* tb, int & b)
01764 {
01765   // Peut-etre rien a faire ?
01766   if (tb->get_etat() == ETATZERO) {
01767     int base_r = b & MSQ_R ;
01768     int base_p = b & MSQ_P ;
01769     switch(base_r){
01770     case(R_CHEBPI_P):
01771       b = R_CHEBPI_I | base_p | T_COSSIN_C ;
01772       break ;
01773     case(R_CHEBPI_I):
01774       b = R_CHEBPI_P | base_p | T_COSSIN_C ;
01775       break ;  
01776     default:
01777       b = base_r | base_p | T_COSSIN_C ;
01778       break;
01779     }
01780     return ;
01781   }
01782   
01783   // Protection
01784   assert(tb->get_etat() == ETATQCQ) ;
01785   
01786   // Pour le confort : nbre de points reels.
01787   int nr = (tb->dim).dim[0] ;
01788   int nt = (tb->dim).dim[1] ;
01789   int np = (tb->dim).dim[2] ;
01790   np = np - 2 ;
01791   
01792   // zone de travail
01793   double* som = new double [nr] ;
01794   
01795   // pt. sur le tableau de double resultat
01796   double* xo = new double[(tb->dim).taille] ;
01797   
01798     // Initialisation a zero :
01799     for (int i=0; i<(tb->dim).taille; i++) {
01800     xo[i] = 0 ; 
01801     }
01802     
01803   // On y va...
01804   double* xi = tb->t ;
01805   double* xci = xi ;    // Pointeurs
01806   double* xco = xo ;    //  courants
01807   
01808   // k = 0
01809   int m = 0 ;       // Parite de l'harmonique en phi
01810   // Dernier j: j = nt-1
01811   // Positionnement
01812   xci += nr * (nt-1) ;
01813   xco += nr * (nt-1) ;
01814     
01815   // Initialisation de som
01816   for (int i=0 ; i<nr ; i++) {
01817     som[i] = 0.5*xci[i] ;
01818     xco[i] = 0. ;
01819   }
01820     
01821   // j suivants
01822   for (int j=nt-2 ; j > 0 ; j--) {
01823     // Positionnement
01824     xci -= 2*nr ;
01825     xco -= nr ;
01826     // Calcul
01827     for (int i=0 ; i<nr ; i++ ) {
01828       som[i] -= 0.5*xci[i] ;
01829       xco[i] = som[i] ;
01830     }   // Fin de la boucle sur r
01831     xci += nr ;
01832     for (int i=0 ; i<nr ; i++ ) {
01833       som[i] = 0.5*xci[i] ;
01834     }
01835   }   // Fin de la boucle sur theta
01836   // j = 0 
01837   xci -= nr ;
01838   xco -= nr ;
01839   for (int i=0; i<nr; i++) {
01840     xco[i] = som[i] ;
01841   }
01842   // Positionnement phi suivant
01843   xci += nr*nt ;
01844   xco += nr*nt ;
01845   
01846   // k=1 
01847   xci += nr*nt ;
01848   xco += nr*nt ;
01849   
01850   for (int k=2 ; k<np+1 ; k++) {
01851     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
01852     
01853     switch(m) {
01854     case 1:     // m impair: cos(impair)
01855       // Dernier j: j = nt-1
01856       // Positionnement
01857       xci += nr * (nt-1) ;
01858       xco += nr * (nt-1) ;
01859       
01860       // Initialisation de som
01861       for (int i=0 ; i<nr ; i++) {
01862     som[i] = -0.5*xci[i] ;
01863     xco[i] = 0.0;
01864       }
01865       
01866       // j suivants
01867       for (int j=nt-2 ; j > 0 ; j--) {
01868     // Positionnement
01869     xci -= 2*nr ;
01870     xco -= nr ;
01871     // Calcul
01872     for (int i=0 ; i<nr ; i++ ) {
01873       som[i] += 0.5*xci[i] ;
01874       xco[i] = som[i] ;
01875     }   // Fin de la boucle sur r
01876     xci +=nr ;
01877     for (int i=0 ; i<nr ; i++ ) {
01878       som[i] = -0.5*xci[i] ;
01879     } 
01880       }   // Fin de la boucle sur theta
01881       xci -= nr ;
01882       for (int i=0; i<nr; i++) {
01883       xco[i] += 0.5*xci[i] ;
01884       }
01885       xco -= nr ;
01886       for (int i=0; i<nr; i++) {
01887     xco[i] = 0. ;
01888       }
01889 
01890       // Positionnement phi suivant
01891       xci += nr*nt ;
01892       xco += nr*nt ;
01893       break ;
01894       
01895     case 0:     // m pair: sin(pair)
01896       // Dernier j: j = nt-1
01897       // Positionnement
01898       xci += nr * (nt-1) ;
01899       xco += nr * (nt-1) ;
01900       
01901       // Initialisation de som
01902       for (int i=0 ; i<nr ; i++) {
01903     som[i] = 0.5*xci[i] ;
01904     xco[i] = 0. ;
01905       }
01906       
01907       // j suivants
01908       for (int j=nt-2 ; j > 0 ; j--) {
01909     // Positionnement
01910     xci -= 2*nr ;
01911     xco -= nr ;
01912     // Calcul
01913     for (int i=0 ; i<nr ; i++ ) {
01914       som[i] -= 0.5*xci[i] ;
01915       xco[i] = som[i] ;
01916     }   // Fin de la boucle sur r
01917     xci += nr ;
01918     for (int i=0 ; i<nr ; i++ ) {
01919       som[i] = 0.5*xci[i] ;
01920     }
01921       }   // Fin de la boucle sur theta
01922       // j = 0 
01923       xci -= nr ;
01924       xco -= nr ;
01925       for (int i=0; i<nr; i++) {
01926     xco[i] = som[i] ;
01927       }
01928       // Positionnement phi suivant
01929       xci += nr*nt ;
01930       xco += nr*nt ;
01931       break ;
01932     }   // Fin du switch
01933   }   // Fin de la boucle en phi
01934   
01935   // On remet les choses la ou il faut
01936   delete [] tb->t ;
01937   tb->t = xo ;
01938   
01939   //Menage
01940   delete [] som ;
01941   
01942   // base de developpement
01943   int base_r = b & MSQ_R ;
01944   int base_p = b & MSQ_P ;
01945   switch(base_r){
01946   case(R_CHEBPI_P):
01947     b = R_CHEBPI_I | base_p | T_COSSIN_C ;
01948     break ;
01949   case(R_CHEBPI_I):
01950     b = R_CHEBPI_P | base_p | T_COSSIN_C ;
01951     break ;  
01952   default:
01953     b = base_r | base_p | T_COSSIN_C ;
01954     break;
01955   }
01956 }
Generated on Tue Feb 7 01:35:18 2012 for LORENE by  doxygen 1.4.6