op_scost.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_scost_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_scost.C,v 1.6 2009/10/10 18:28:11 j_novak Exp $" ;
00024 
00025 /*
00026  * $Id: op_scost.C,v 1.6 2009/10/10 18:28:11 j_novak Exp $
00027  * $Log: op_scost.C,v $
00028  * Revision 1.6  2009/10/10 18:28:11  j_novak
00029  * New bases T_COS and T_SIN.
00030  *
00031  * Revision 1.5  2007/12/21 10:43:38  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:40  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:36  eric
00050  * Initialisation a zero du tableau xo avant le calcul.
00051  *
00052  * Revision 1.2  1999/11/22  14:34:16  novak
00053  * Suppression des variables inutiles
00054  *
00055  * Revision 1.1  1999/11/16  13:38:00  novak
00056  * Initial revision
00057  *
00058  *
00059  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_scost.C,v 1.6 2009/10/10 18:28:11 j_novak Exp $
00060  *
00061  */
00062 
00063 /* 
00064  * Ensemble des routines de base de division par rapport a cost theta
00065  * (Utilisation interne)
00066  * 
00067  *  void _scost_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 _scost_pas_prevu(Tbl * tb, int& base) {
00083     cout << "scost 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 _scost_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     // pt. sur le tableau de double resultat
00124     double* somP = new double [nr] ;
00125     double* somN = new double [nr] ;
00126     
00127     // pt. sur le tableau de double resultat
00128     double* xo = new double[(tb->dim).taille] ;
00129     
00130     // Initialisation a zero :
00131     for (int i=0; i<(tb->dim).taille; i++) {
00132     xo[i] = 0 ; 
00133     }
00134     
00135     // On y va...
00136     double* xi = tb->t ;
00137     double* xci = xi ;  // Pointeurs
00138     double* xco = xo ;  //  courants
00139     
00140     // k = 0
00141         
00142     // Dernier j: j = nt-1
00143     // Positionnement
00144     xci += nr * (nt-1) ;
00145     xco += nr * (nt-1) ;
00146     
00147     // Initialisation de som 
00148     for (int i=0 ; i<nr ; i++) {
00149         somP[i] = 0. ;
00150         somN[i] = 0. ;
00151         xco[i] = 0. ;   // mise a zero du dernier coef en theta
00152     }
00153     
00154     // j suivants
00155     for (int j=nt-2 ; j >= 0 ; j--) {
00156         // Positionnement
00157         xci -= nr ;
00158         xco -= nr ;
00159         // Calcul
00160         for (int i=0 ; i<nr ; i++ ) {
00161           if((j%2) == 1) {
00162         somN[i] = -somN[i] ;
00163         somN[i] += 2*xci[i] ;
00164         xco[i] = somP[i] ;
00165           }
00166           else {
00167         somP[i] = -somP[i] ;
00168         somP[i] += 2*xci[i] ; 
00169         xco[i] = somN[i] ;
00170           }
00171         }   // Fin de la boucle sur r
00172     }   // Fin de la boucle sur theta
00173     // j=0 : le facteur 2...
00174     for (int i=0 ; i<nr ; i++) xco[i] *= .5 ;
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         somP[i] = 0. ;
00195         somN[i] = 0. ;
00196         xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
00197     }
00198     
00199     // j suivants
00200     for (int j=nt-2 ; j >= 0 ; j--) {
00201         // Positionnement
00202         xci -= nr ;
00203         xco -= nr ;
00204         // Calcul
00205         for (int i=0 ; i<nr ; i++ ) {
00206         if((j%2) == 1) {
00207             somN[i] = -somN[i];
00208             somN[i] += 2 * xci[i] ;
00209             xco[i] = somP[i] ;
00210         }
00211         else {
00212             somP[i] = -somP[i];
00213             somP[i] += 2 * xci[i] ;
00214             xco[i] = somN[i];
00215         }
00216         }   // Fin de la boucle sur r
00217     }   // Fin de la boucle sur theta
00218     // j=0 : facteur 2 ...
00219     for (int i=0 ; i<nr ; i++) xco[i] *= 0.5 ;
00220 
00221     // Positionnement phi suivant
00222     xci += nr*nt ;
00223     xco += nr*nt ;
00224     }   // Fin de la boucle sur phi
00225 
00226     // On remet les choses la ou il faut
00227     delete [] tb->t ;
00228     tb->t = xo ;
00229   
00230     //Menage
00231     delete [] somP ;
00232     delete [] somN ;
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_COS ;
00240         break ;
00241         case(R_CHEBPI_I):
00242         b = R_CHEBPI_P | base_p | T_COS ;
00243         break ;  
00244         default:
00245         b = base_r | base_p | T_COS ;
00246         break;
00247     }
00248 }
00249             
00250             //------------
00251             // cas T_SIN
00252             //------------
00253             
00254 void _scost_t_sin(Tbl* tb, int & b)
00255 {
00256 
00257 
00258     // Peut-etre rien a faire ?
00259     if (tb->get_etat() == ETATZERO) {
00260     int base_r = b & MSQ_R ;
00261     int base_p = b & MSQ_P ;
00262     switch(base_r){
00263         case(R_CHEBPI_P):
00264         b = R_CHEBPI_I | base_p | T_SIN ;
00265         break ;
00266         case(R_CHEBPI_I):
00267         b = R_CHEBPI_P | base_p | T_SIN ;
00268         break ;  
00269         default:
00270         b = base_r | base_p | T_SIN ;
00271         break;
00272     }
00273     return ;
00274     }
00275     
00276     // Protection
00277     assert(tb->get_etat() == ETATQCQ) ;
00278     
00279     // Pour le confort : nbre de points reels.
00280     int nr = (tb->dim).dim[0] ;
00281     int nt = (tb->dim).dim[1] ;
00282     int np = (tb->dim).dim[2] ;
00283     np = np - 2 ;
00284     
00285     // zone de travail
00286     double* somP = new double [nr] ;
00287     double* somN = new double [nr] ;
00288     
00289     // pt. sur le tableau de double resultat
00290     double* xo = new double[(tb->dim).taille] ;
00291     
00292     // Initialisation a zero :
00293     for (int i=0; i<(tb->dim).taille; i++) {
00294     xo[i] = 0 ; 
00295     }
00296     
00297     // On y va...
00298     double* xi = tb->t ;
00299     double* xci = xi ;  // Pointeurs
00300     double* xco = xo ;  //  courants
00301     
00302     // k = 0
00303     
00304     // Dernier j: j = nt-1
00305     // Positionnement
00306     xci += nr * (nt-1) ;
00307     xco += nr * (nt-1) ;
00308     
00309     // Initialisation de som
00310     for (int i=0 ; i<nr ; i++) {
00311         somP[i] = 0. ;
00312         somN[i] = 0. ;
00313         xco[i] = 0. ;
00314     }
00315     
00316     // j suivants
00317     for (int j=nt-2 ; j >= 0 ; j--) {
00318       // Positionnement
00319       xci -= nr ;
00320       xco -= nr ;
00321       // Calcul
00322       for (int i=0 ; i<nr ; i++ ) {
00323     if((j%2) == 1) {
00324       somN[i] = -somN[i] ;
00325       somN[i] += 2*xci[i] ;
00326       xco[i] = somP[i] ;
00327     }
00328     else {
00329       somP[i] = -somP[i] ;
00330       somP[i] += 2*xci[i] ;       
00331       xco[i] = somN[i] ;
00332     }
00333       } // Fin de la boucle sur r
00334     }   // Fin de la boucle sur theta
00335     // j=0 : sin(0*theta)
00336     for (int i=0 ; i<nr ; i++) xco[i] = 0. ;
00337 
00338     // Positionnement phi suivant
00339     xci += nr*nt ;
00340     xco += nr*nt ;
00341 
00342     // k = 1
00343     xci += nr*nt ;
00344     xco += nr*nt ;
00345     
00346     // k >= 2
00347     for (int k=2 ; k<np+1 ; k++) {
00348     
00349       // Dernier j: j = nt-1
00350       // Positionnement
00351       xci += nr * (nt-1) ;
00352       xco += nr * (nt-1) ;
00353     
00354     // Initialisation de som
00355     for (int i=0 ; i<nr ; i++) {
00356         somP[i] = 0. ;
00357         somN[i] = 0. ;
00358         xco[i] = 0. ;
00359     }
00360     
00361     // j suivants
00362     for (int j=nt-2 ; j >= 0 ; j--) {
00363         // Positionnement
00364         xci -= nr ;
00365         xco -= nr ;
00366         // Calcul
00367         for (int i=0 ; i<nr ; i++ ) {
00368         if((j%2) == 1) {
00369           somN[i] = -somN[i] ;
00370           somN[i] += 2*xci[i] ;
00371           xco[i] = somP[i] ;
00372           }
00373           else {
00374           somP[i] = -somP[i] ;
00375           somP[i] += 2*xci[i] ;       
00376           xco[i] = somN[i] ;
00377           }
00378         }   // Fin de la boucle sur r
00379     }   // Fin de la boucle sur theta
00380     // j=0 : sin(0*theta)
00381     for (int i=0 ; i<nr ; i++) xco[i] = 0. ;
00382     
00383     // Positionnement phi suivant
00384     xci += nr*nt ;
00385     xco += nr*nt ;
00386 
00387     }   // Fin de la boucle sur phi
00388 
00389     // On remet les choses la ou il faut
00390     delete [] tb->t ;
00391     tb->t = xo ;
00392     
00393     //Menage
00394     delete [] somP ;
00395     delete [] somN ;
00396 
00397     // base de developpement
00398     int base_r = b & MSQ_R ;
00399     int base_p = b & MSQ_P ;
00400     switch(base_r){
00401     case(R_CHEBPI_P):
00402         b = R_CHEBPI_I | base_p | T_SIN ;
00403         break ;
00404     case(R_CHEBPI_I):
00405         b = R_CHEBPI_P | base_p | T_SIN ;
00406         break ;  
00407     default:
00408         b = base_r | base_p | T_SIN ;
00409         break;
00410     }
00411 }
00412             //----------------
00413             // cas T_COS_P
00414             //----------------
00415             
00416 void _scost_t_cos_p(Tbl* tb, int & b)
00417 {
00418 
00419     // Peut-etre rien a faire ?
00420     if (tb->get_etat() == ETATZERO) {
00421     int base_r = b & MSQ_R ;
00422     int base_p = b & MSQ_P ;
00423     b = base_r | base_p | T_COS_I ;
00424     return ;
00425     }
00426     
00427     // Protection
00428     assert(tb->get_etat() == ETATQCQ) ;
00429     
00430     // Pour le confort : nbre de points reels.
00431     int nr = (tb->dim).dim[0] ;
00432     int nt = (tb->dim).dim[1] ;
00433     int np = (tb->dim).dim[2] ;
00434     np = np - 2 ;
00435     
00436     // zone de travail
00437     double* som = new double [nr] ;
00438     
00439     // pt. sur le tableau de double resultat
00440     double* xo = new double[(tb->dim).taille] ;
00441     
00442     // Initialisation a zero :
00443     for (int i=0; i<(tb->dim).taille; i++) {
00444     xo[i] = 0 ; 
00445     }
00446     
00447     // On y va...
00448     double* xi = tb->t ;
00449     double* xci = xi ;  // Pointeurs
00450     double* xco = xo ;  //  courants
00451     
00452     // k = 0
00453         
00454     // Dernier j: j = nt-1
00455     // Positionnement
00456     xci += nr * (nt) ;
00457     xco += nr * (nt-1) ;
00458     
00459     // Initialisation de som 
00460     for (int i=0 ; i<nr ; i++) {
00461         som[i] = 0. ;
00462         xco[i] = 0. ;   // mise a zero du dernier coef en theta
00463     }
00464     
00465     // j suivants
00466     for (int j=nt-2 ; j >= 0 ; j--) {
00467         // Positionnement
00468         xci -= nr ;
00469         xco -= nr ;
00470         // Calcul
00471         for (int i=0 ; i<nr ; i++ ) {
00472         som[i] += 2*xci[i] ;
00473         xco[i] = som[i] ;
00474         som[i] = -som[i] ;
00475         }   // Fin de la boucle sur r
00476     }   // Fin de la boucle sur theta
00477     // Positionnement phi suivant
00478     xci -= nr ;
00479     xci += nr*nt ;
00480     xco += nr*nt ;
00481 
00482     // k = 1
00483     xci += nr*nt ;
00484     xco += nr*nt ;
00485     
00486     // k >= 2
00487     for (int k=2 ; k<np+1 ; k++) {
00488     
00489     // Dernier j: j = nt-1
00490     // Positionnement
00491     xci += nr * (nt) ;
00492     xco += nr * (nt-1) ;
00493     
00494     // Initialisation de som
00495     for (int i=0 ; i<nr ; i++) {
00496         som[i] = 0. ;
00497         xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
00498     }
00499     
00500     // j suivants
00501     for (int j=nt-2 ; j >= 0 ; j--) {
00502         // Positionnement
00503         xci -= nr ;
00504         xco -= nr ;
00505         // Calcul
00506         for (int i=0 ; i<nr ; i++ ) {
00507         som[i] += 2*xci[i] ;
00508         xco[i] = som[i] ;
00509         som[i] = - som[i] ;
00510         }   // Fin de la boucle sur r
00511     }   // Fin de la boucle sur theta
00512     // Positionnement phi suivant
00513     xci -= nr ;
00514     xci += nr*nt ;
00515     xco += nr*nt ;
00516     }   // Fin de la boucle sur phi
00517 
00518     // On remet les choses la ou il faut
00519     delete [] tb->t ;
00520     tb->t = xo ;
00521   
00522     //Menage
00523     delete [] som ;
00524 
00525     // base de developpement
00526     int base_r = b & MSQ_R ;
00527     int base_p = b & MSQ_P ;
00528     b = base_r | base_p | T_COS_I ;
00529 }
00530             
00531             //--------------
00532             // cas T_SIN_P
00533             //--------------
00534             
00535 void _scost_t_sin_p(Tbl* tb, int & b)
00536 {
00537 
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 ;
00577     xco += nr * (nt-1) ;
00578     
00579     // Initialisation de som
00580     for (int i=0 ; i<nr ; i++) {
00581         som[i] = 0. ;
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] += 2*xci[i] ;
00593         xco[i] = som[i] ;
00594         som[i] = -som[i] ;
00595         }   // Fin de la boucle sur r
00596     }   // Fin de la boucle sur theta
00597     // Positionnement phi suivant
00598     xci -= nr ;
00599     xci += nr*nt ;
00600     xco += nr*nt ;
00601 
00602     // k = 1
00603     xci += nr*nt ;
00604     xco += nr*nt ;
00605     
00606     // k >= 2
00607     for (int k=2 ; k<np+1 ; k++) {
00608     
00609     // Dernier j: j = nt-1
00610     // Positionnement
00611     xci += nr * nt ;
00612     xco += nr * (nt-1) ;
00613     
00614     // Initialisation de som
00615     for (int i=0 ; i<nr ; i++) {
00616         som[i] = 0. ;
00617         xco[i] = 0. ;
00618     }
00619     
00620     // j suivants
00621     for (int j=nt-2 ; j >= 0 ; j--) {
00622         // Positionnement
00623         xci -= nr ;
00624         xco -= nr ;
00625         // Calcul
00626         for (int i=0 ; i<nr ; i++ ) {
00627         som[i] += 2*xci[i] ;
00628         xco[i] = som[i] ;
00629         som[i] = -som[i] ;
00630         }   // Fin de la boucle sur r
00631     }   // Fin de la boucle sur theta
00632     // Positionnement phi suivant
00633     xci -= nr ;
00634     xci += nr*nt ;
00635     xco += nr*nt ;
00636     }   // Fin de la boucle sur phi
00637 
00638     // On remet les choses la ou il faut
00639     delete [] tb->t ;
00640     tb->t = xo ;
00641     
00642     //Menage
00643     delete [] som ;
00644 
00645     // base de developpement
00646     int base_r = b & MSQ_R ;
00647     int base_p = b & MSQ_P ;
00648     b = base_r | base_p | T_SIN_I ;
00649    
00650 }
00651             //--------------
00652             // cas T_SIN_I
00653             //--------------
00654             
00655 void _scost_t_sin_i(Tbl* tb, int & b)
00656 {
00657 
00658 
00659     // Peut-etre rien a faire ?
00660     if (tb->get_etat() == ETATZERO) {
00661     int base_r = b & MSQ_R ;
00662     int base_p = b & MSQ_P ;
00663     b = base_r | base_p | T_SIN_P ;
00664     return ;
00665     }
00666     
00667     // Protection
00668     assert(tb->get_etat() == ETATQCQ) ;
00669     
00670     // Pour le confort : nbre de points reels.
00671     int nr = (tb->dim).dim[0] ;
00672     int nt = (tb->dim).dim[1] ;
00673     int np = (tb->dim).dim[2] ;
00674     np = np - 2 ;
00675     
00676     // zone de travail
00677     double* som = new double [nr] ;
00678     
00679     // pt. sur le tableau de double resultat
00680     double* xo = new double[(tb->dim).taille] ;
00681     
00682     // Initialisation a zero :
00683     for (int i=0; i<(tb->dim).taille; i++) {
00684     xo[i] = 0 ; 
00685     }
00686     
00687     // On y va...
00688     double* xi = tb->t ;
00689     double* xci = xi ;  // Pointeurs
00690     double* xco = xo ;  //  courants
00691     
00692        
00693     // k = 0
00694     
00695     // Dernier j: j = nt-1
00696     // Positionnement
00697     xci += nr * (nt-1) ;
00698     xco += nr * (nt-1) ;
00699 
00700     // Initialisation de som
00701     for (int i=0 ; i<nr ; i++) {
00702         xco[i] = 0. ;
00703         som[i] = 0. ;
00704     }
00705     
00706     // j suivants
00707     for (int j=nt-2 ; j >= 1 ; j--) {
00708         // Positionnement
00709         xci -= nr ;
00710         xco -= nr ;
00711         // Calcul
00712         for (int i=0 ; i<nr ; i++ ) {
00713         som[i] += 2*xci[i] ;
00714         xco[i] = som[i] ;
00715         som[i] = -som[i] ;
00716         }   // Fin de la boucle sur r
00717     }   // Fin de la boucle sur theta
00718     if (nt > 1) {
00719         xci -= nr ;
00720         xco -= nr ;
00721         // Calcul pour le premier theta
00722         for (int i=0 ; i<nr ; i++ ) {
00723         xco[i] =0. ;
00724         }
00725     }
00726     // Positionnement phi suivant
00727     xci += nr*nt ;
00728     xco += nr*nt ;
00729 
00730     // k = 1
00731     xci += nr*nt ;
00732     xco += nr*nt ;
00733     
00734     // k >= 2
00735     for (int k=2 ; k<np+1 ; k++) {
00736     
00737     // Dernier j: j = nt-1
00738     // Positionnement
00739     xci += nr * (nt-1) ;
00740     xco += nr * (nt-1) ;
00741 
00742     // Initialisation de som
00743     for (int i=0 ; i<nr ; i++) {
00744         xco[i] = 0. ;
00745         som[i] = 0. ;
00746     }
00747     
00748     // j suivants
00749     for (int j=nt-2 ; j >= 1 ; j--) {
00750         // Positionnement
00751         xci -= nr ;
00752         xco -= nr ;
00753         // Calcul
00754         for (int i=0 ; i<nr ; i++ ) {
00755         som[i] += 2*xci[i] ;
00756         xco[i] = som[i] ;
00757         som[i] = -som[i] ;
00758         }   // Fin de la boucle sur r
00759     }   // Fin de la boucle sur theta
00760     xci -= nr ;
00761     xco -= nr ;
00762     // Calcul pour le premier theta
00763     for (int i=0 ; i<nr ; i++ ) {
00764       xco[i] =0. ;
00765     }
00766     // Positionnement phi suivant
00767     xci += nr*nt ;
00768     xco += nr*nt ;
00769     }   // Fin de la boucle sur phi
00770 
00771 
00772     // On remet les choses la ou il faut
00773     delete [] tb->t ;
00774     tb->t = xo ;
00775     
00776     //Menage
00777     delete [] som ;
00778 
00779     // base de developpement
00780     int base_r = b & MSQ_R ;
00781     int base_p = b & MSQ_P ;
00782     b = base_r | base_p | T_SIN_P ;
00783 
00784 }
00785             //---------------------
00786             // cas T_COS_I
00787             //---------------------
00788 void _scost_t_cos_i(Tbl* tb, int & b)
00789 {
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         xco[i] = 0. ;   // mise a zero du dernier coef en theta
00835     }
00836     
00837     // j suivants
00838     for (int j=nt-2 ; j >= 0 ; j--) {
00839         // Positionnement
00840         xci -= nr ;
00841         xco -= nr ;
00842         // Calcul
00843         for (int i=0 ; i<nr ; i++ ) {
00844         som[i] += 2*xci[i] ;
00845         xco[i] = som[i] ;
00846         som[i] = - som[i] ;
00847         }   // Fin de la boucle sur r
00848     }   // Fin de la boucle sur theta
00849     // Normalisation du premier theta
00850     for (int i=0 ; i<nr ; i++) {
00851         xco[i] *= .5 ;
00852     }
00853     // Positionnement phi suivant
00854     xci += nr*nt ;
00855     xco += nr*nt ;
00856 
00857     // k = 1
00858     xci += nr*nt ;
00859     xco += nr*nt ;
00860     
00861     // k >= 2
00862     for (int k=2 ; k<np+1 ; k++) {
00863     
00864     // Dernier j: j = nt-1
00865     // Positionnement
00866     xci += nr * (nt-1) ;
00867     xco += nr * (nt-1) ;
00868     
00869     // Initialisation de som
00870     for (int i=0 ; i<nr ; i++) {
00871         som[i] = 0. ;
00872         xco[i] = 0. ;   // mise a zero du dernier coef en theta
00873     }
00874     
00875     // j suivants
00876     for (int j=nt-2 ; j >= 0 ; j--) {
00877         // Positionnement
00878         xci -= nr ;
00879         xco -= nr ;
00880         // Calcul
00881         for (int i=0 ; i<nr ; i++ ) {
00882         som[i] += 2*xci[i] ;
00883         xco[i] = som[i] ;
00884         som[i] = -som[i] ;
00885         }   // Fin de la boucle sur r
00886     }   // Fin de la boucle sur theta
00887     // Normalisation du premier theta
00888     for (int i=0 ; i<nr ; i++) {
00889         xco[i] *= .5 ;
00890     }
00891     // Positionnement phi suivant
00892     xci += nr*nt ;
00893     xco += nr*nt ;
00894     }   // Fin de la boucle sur phi
00895 
00896     // On remet les choses la ou il faut
00897     delete [] tb->t ;
00898     tb->t = xo ;
00899     
00900     //Menage
00901     delete [] som ;
00902 
00903     // base de developpement
00904     int base_r = b & MSQ_R ;
00905     int base_p = b & MSQ_P ;
00906     b = base_r | base_p | T_COS_P ;
00907   
00908 }
00909             //----------------------
00910             // cas T_COSSIN_CP
00911             //----------------------
00912 void _scost_t_cossin_cp(Tbl* tb, int & b)
00913 {
00914 
00915     // Peut-etre rien a faire ?
00916     if (tb->get_etat() == ETATZERO) {
00917     int base_r = b & MSQ_R ;
00918     int base_p = b & MSQ_P ;
00919     b = base_r | base_p | T_COSSIN_CI ;
00920     return ;
00921     }
00922     
00923     // Protection
00924     assert(tb->get_etat() == ETATQCQ) ;
00925     
00926     // Pour le confort : nbre de points reels.
00927     int nr = (tb->dim).dim[0] ;
00928     int nt = (tb->dim).dim[1] ;
00929     int np = (tb->dim).dim[2] ;
00930     np = np - 2 ;
00931     
00932     // zone de travail
00933     double* som = new double [nr] ;
00934     
00935     // pt. sur le tableau de double resultat
00936     double* xo = new double[(tb->dim).taille] ;
00937     
00938     // Initialisation a zero :
00939     for (int i=0; i<(tb->dim).taille; i++) {
00940     xo[i] = 0 ; 
00941     }
00942     
00943     // On y va...
00944     double* xi = tb->t ;
00945     double* xci = xi ;  // Pointeurs
00946     double* xco = xo ;  //  courants
00947     
00948     // k = 0
00949     int m = 0 ;     // Parite de l'harmonique en phi
00950     
00951     // Dernier j: j = nt-1
00952     // Positionnement
00953     xci += nr * (nt) ;
00954     xco += nr * (nt-1) ;
00955     
00956     // Initialisation de som 
00957     for (int i=0 ; i<nr ; i++) {
00958         som[i] = 0. ;
00959         xco[i] = 0. ;   // mise a zero du dernier coef en theta
00960     }
00961     
00962     // j suivants
00963     for (int j=nt-2 ; j >= 0 ; j--) {
00964         // Positionnement
00965         xci -= nr ;
00966         xco -= nr ;
00967         // Calcul
00968         for (int i=0 ; i<nr ; i++ ) {
00969         som[i] += 2*xci[i] ;
00970         xco[i] = som[i] ;
00971         som[i] = -som[i] ;
00972         }   // Fin de la boucle sur r
00973     }   // Fin de la boucle sur theta
00974     // Positionnement phi suivant
00975     xci -= nr ;
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) ;
00992     xco += nr * (nt-1) ;
00993     
00994     // Initialisation de som 
00995     for (int i=0 ; i<nr ; i++) {
00996       som[i] = 0. ;
00997       xco[i] = 0. ; // mise a zero du dernier coef 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] += 2*xci[i] ;
01008         xco[i] = som[i] ;
01009         som[i] = -som[i] ;
01010       } // Fin de la boucle sur r
01011     }   // Fin de la boucle sur theta
01012     // Positionnement phi suivant
01013     xci -= nr ;
01014     xci += nr*nt ;
01015     xco += nr*nt ;
01016     break ;
01017           
01018       case 1:       // m impair: sin(impair)
01019     // Dernier j: j = nt-1
01020     // Positionnement
01021     xci += nr * (nt-1) ;
01022     xco += nr * (nt-1) ;
01023 
01024     // Initialisation de som
01025     for (int i=0 ; i<nr ; i++) {
01026       xco[i] = 0. ;
01027       som[i] = 0. ;
01028     }
01029     
01030     // j suivants
01031     for (int j=nt-2 ; j >= 1 ; j--) {
01032       // Positionnement
01033       xci -= nr ;
01034       xco -= nr ;
01035       // Calcul
01036       for (int i=0 ; i<nr ; i++ ) {
01037         som[i] += 2*xci[i] ;
01038         xco[i] = som[i] ;
01039         som[i] = -som[i] ;
01040       } // Fin de la boucle sur r
01041     }   // Fin de la boucle sur theta
01042     xci -= nr ;
01043     xco -= nr ;
01044     // Calcul pour le premier theta
01045     for (int i=0 ; i<nr ; i++ ) {
01046       xco[i] =0. ;
01047     }
01048     // Positionnement phi suivant
01049     xci += nr*nt ;
01050     xco += nr*nt ;
01051     break;
01052       } // Fin du switch
01053     }   // Fin de la boucle sur phi
01054 
01055     // On remet les choses la ou il faut
01056     delete [] tb->t ;
01057     tb->t = xo ;
01058     
01059     //Menage
01060     delete [] som ;
01061 
01062     // base de developpement
01063     int base_r = b & MSQ_R ;
01064     int base_p = b & MSQ_P ;
01065     b = base_r | base_p | T_COSSIN_CI ;
01066 }
01067             
01068             //------------------
01069             // cas T_COSSIN_CI
01070             //----------------
01071 void _scost_t_cossin_ci(Tbl* tb, int & b)
01072 {
01073     // Peut-etre rien a faire ?
01074     if (tb->get_etat() == ETATZERO) {
01075     int base_r = b & MSQ_R ;
01076     int base_p = b & MSQ_P ;
01077     b = base_r | base_p | T_COSSIN_CP ;
01078     return ;
01079     }
01080     
01081     // Protection
01082     assert(tb->get_etat() == ETATQCQ) ;
01083     
01084     // Pour le confort : nbre de points reels.
01085     int nr = (tb->dim).dim[0] ;
01086     int nt = (tb->dim).dim[1] ;
01087     int np = (tb->dim).dim[2] ;
01088     np = np - 2 ;
01089     
01090     // zone de travail
01091     double* som = new double [nr] ;
01092     
01093     // pt. sur le tableau de double resultat
01094     double* xo = new double[(tb->dim).taille] ;
01095     
01096     // Initialisation a zero :
01097     for (int i=0; i<(tb->dim).taille; i++) {
01098     xo[i] = 0 ; 
01099     }
01100     
01101     // On y va...
01102     double* xi = tb->t ;
01103     double* xci = xi ;  // Pointeurs
01104     double* xco = xo ;  //  courants
01105 
01106     // k = 0
01107     int m = 0 ;     // Parite de l'harmonique en phi
01108     
01109     
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         xco[i] = 0. ;   // mise a zero du dernier coef en theta
01119     }
01120     
01121     // j suivants
01122     for (int j=nt-2 ; j >= 0 ; j--) {
01123         // Positionnement
01124         xci -= nr ;
01125         xco -= nr ;
01126         // Calcul
01127         for (int i=0 ; i<nr ; i++ ) {
01128         som[i] += 2*xci[i] ;
01129         xco[i] = som[i] ;
01130         som[i] = - som[i] ;
01131         }   // Fin de la boucle sur r
01132     }   // Fin de la boucle sur theta
01133     // Normalisation du premier theta
01134     for (int i=0 ; i<nr ; i++) {
01135         xco[i] *= .5 ;
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       xco[i] = 0. ; // mise a zero du dernier coef en theta
01160     }
01161           
01162     // j suivants
01163     for (int j=nt-2 ; j >= 0 ; j--) {
01164       // Positionnement
01165       xci -= nr ;
01166       xco -= nr ;
01167       // Calcul
01168       for (int i=0 ; i<nr ; i++ ) {
01169         som[i] += 2*xci[i] ;
01170         xco[i] = som[i] ;
01171         som[i] = - som[i] ;
01172       } // Fin de la boucle sur r
01173     }   // Fin de la boucle sur theta
01174     // Normalisation du premier theta
01175     for (int i=0 ; i<nr ; i++) {
01176       xco[i] *= .5 ;
01177     }
01178     // Positionnement phi suivant
01179     xci += nr*nt ;
01180     xco += nr*nt ;
01181     break ;
01182     
01183       case 1:
01184     // m impair: sin(pair)
01185     // Dernier j: j = nt-1
01186     // Positionnement
01187     xci += nr * nt ;
01188     xco += nr * (nt-1) ;
01189     
01190     // Initialisation de som
01191     for (int i=0 ; i<nr ; i++) {
01192       som[i] = 0. ;
01193       xco[i] = 0. ;
01194     }
01195     
01196     // j suivants
01197     for (int j=nt-2 ; j >= 0 ; j--) {
01198       // Positionnement
01199       xci -= nr ;
01200       xco -= nr ;
01201       // Calcul
01202       for (int i=0 ; i<nr ; i++ ) {
01203         som[i] += 2*xci[i] ;
01204         xco[i] = som[i] ;
01205         som[i] = -som[i] ;
01206       } // Fin de la boucle sur r
01207     }   // Fin de la boucle sur theta
01208     // Positionnement phi suivant
01209     xci -= nr ;
01210     xci += nr*nt ;
01211     xco += nr*nt ;
01212     break ;
01213       }   // Fin du switch
01214     }   // Fin de la boucle en phi
01215     
01216     // On remet les choses la ou il faut
01217     delete [] tb->t ;
01218     tb->t = xo ;
01219     
01220     //Menage
01221     delete [] som ;
01222 
01223     // base de developpement
01224     int base_r = b & MSQ_R ;
01225     int base_p = b & MSQ_P ;
01226     b = base_r | base_p | T_COSSIN_CP ;
01227 }
01228 
01229             //--------------------- 
01230             // cas T_COSSIN_SI
01231             //----------------
01232 void _scost_t_cossin_si(Tbl* tb, int & b)
01233 {
01234   // Peut-etre rien a faire ?
01235   if (tb->get_etat() == ETATZERO) {
01236     int base_r = b & MSQ_R ;
01237     int base_p = b & MSQ_P ;
01238     b = base_r | base_p | T_COSSIN_SP ;
01239     return ;
01240   }
01241     
01242   // Protection
01243   assert(tb->get_etat() == ETATQCQ) ;
01244   
01245   // Pour le confort : nbre de points reels.
01246   int nr = (tb->dim).dim[0] ;
01247   int nt = (tb->dim).dim[1] ;
01248   int np = (tb->dim).dim[2] ;
01249   np = np - 2 ;
01250   
01251   // zone de travail
01252   double* som = new double [nr] ;
01253   
01254   // pt. sur le tableau de double resultat
01255   double* xo = new double[(tb->dim).taille] ;
01256   
01257     // Initialisation a zero :
01258     for (int i=0; i<(tb->dim).taille; i++) {
01259     xo[i] = 0 ; 
01260     }
01261     
01262   // On y va...
01263   double* xi = tb->t ;
01264   double* xci = xi ;    // Pointeurs
01265   double* xco = xo ;    //  courants
01266   
01267   // k = 0
01268   int m = 0 ;       // Parite de l'harmonique en phi
01269   
01270   // Dernier j: j = nt-1
01271   // Positionnement
01272   xci += nr * (nt-1) ;
01273   xco += nr * (nt-1) ;
01274   
01275   // Initialisation de som
01276   for (int i=0 ; i<nr ; i++) {
01277     xco[i] = 0. ;
01278     som[i] = 0. ;
01279   }
01280   
01281   // j suivants
01282   for (int j=nt-2 ; j >= 1 ; j--) {
01283     // Positionnement
01284     xci -= nr ;
01285     xco -= nr ;
01286     // Calcul
01287     for (int i=0 ; i<nr ; i++ ) {
01288       som[i] += 2*xci[i] ;
01289       xco[i] = som[i] ;
01290       som[i] = -som[i] ;
01291     }   // Fin de la boucle sur r
01292   }   // Fin de la boucle sur theta
01293   xci -= nr ;
01294   xco -= nr ;
01295   // Calcul pour le premier theta
01296   for (int i=0 ; i<nr ; i++ ) {
01297     xco[i] =0. ;
01298   }
01299   // Positionnement phi suivant
01300   xci += nr*nt ;
01301   xco += nr*nt ;
01302   
01303   // k=1 
01304   xci += nr*nt ;
01305   xco += nr*nt ;
01306   
01307   // k >= 2
01308   for (int k=2 ; k<np+1 ; k++) {
01309     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
01310     
01311     switch(m) {
01312     case 0:     // m pair: sin(impair)
01313       // Dernier j: j = nt-1
01314       // Positionnement
01315       xci += nr * (nt-1) ;
01316       xco += nr * (nt-1) ;
01317       
01318       // Initialisation de som
01319       for (int i=0 ; i<nr ; i++) {
01320     xco[i] = 0. ;
01321     som[i] = 0. ;
01322       }
01323       
01324       // j suivants
01325       for (int j=nt-2 ; j >= 1 ; j--) {
01326     // Positionnement
01327     xci -= nr ;
01328     xco -= nr ;
01329     // Calcul
01330     for (int i=0 ; i<nr ; i++ ) {
01331       som[i] += 2*xci[i] ;
01332       xco[i] = som[i] ;
01333       som[i] = -som[i] ;
01334     }   // Fin de la boucle sur r
01335       }   // Fin de la boucle sur theta
01336       xci -= nr ;
01337       xco -= nr ;
01338       // Calcul pour le premier theta
01339       for (int i=0 ; i<nr ; i++ ) {
01340     xco[i] =0. ;
01341       }
01342       // Positionnement phi suivant
01343       xci += nr*nt ;
01344       xco += nr*nt ;
01345       break ;
01346       
01347     case 1: // m impair cos(pair)
01348       // Dernier j: j = nt-1
01349       // Positionnement
01350       xci += nr * (nt) ;
01351       xco += nr * (nt-1) ;
01352       
01353       // Initialisation de som
01354       for (int i=0 ; i<nr ; i++) {
01355     som[i] = 0. ;
01356     xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
01357       }
01358       
01359       // j suivants
01360       for (int j=nt-2 ; j >= 0 ; j--) {
01361     // Positionnement
01362     xci -= nr ;
01363     xco -= nr ;
01364     // Calcul
01365     for (int i=0 ; i<nr ; i++ ) {
01366       som[i] += 2*xci[i] ;
01367       xco[i] = som[i] ;
01368       som[i] = - som[i] ;
01369     }   // Fin de la boucle sur r
01370       }   // Fin de la boucle sur theta
01371       // Positionnement phi suivant
01372       xci -= nr ;
01373       xci += nr*nt ;
01374       xco += nr*nt ;
01375       
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 _scost_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   
01435   // Dernier j: j = nt-1
01436   // Positionnement
01437   xci += nr * nt ;
01438   xco += nr * (nt-1) ;
01439   
01440   // Initialisation de som
01441   for (int i=0 ; i<nr ; i++) {
01442     som[i] = 0. ;
01443     xco[i] = 0. ;
01444   }
01445     
01446   // j suivants
01447   for (int j=nt-2 ; j >= 0 ; j--) {
01448     // Positionnement
01449     xci -= nr ;
01450     xco -= nr ;
01451     // Calcul
01452     for (int i=0 ; i<nr ; i++ ) {
01453       som[i] += 2*xci[i] ;
01454       xco[i] = som[i] ;
01455       som[i] = -som[i] ;
01456     }   // Fin de la boucle sur r
01457   }   // Fin de la boucle sur theta
01458   // Positionnement phi suivant
01459   xci -= nr ;
01460   xci += nr*nt ;
01461   xco += nr*nt ;
01462   
01463   // k=1 
01464   xci += nr*nt ;
01465   xco += nr*nt ;
01466   
01467   for (int k=2 ; k<np+1 ; k++) {
01468     m = (k/2) % 2 ;     // Parite de l'harmonique en phi
01469     
01470     switch(m) {
01471     case 1:     // m impair: cos(impair)
01472       // Dernier j: j = nt-1
01473       // Positionnement
01474       xci += nr * (nt-1) ;
01475       xco += nr * (nt-1) ;
01476       
01477       // Initialisation de som
01478       for (int i=0 ; i<nr ; i++) {
01479     som[i] = 0. ;
01480     xco[i] = 0. ;   // mise a zero du dernier coef en theta
01481       }
01482       
01483       // j suivants
01484       for (int j=nt-2 ; j >= 0 ; j--) {
01485     // Positionnement
01486     xci -= nr ;
01487     xco -= nr ;
01488     // Calcul
01489     for (int i=0 ; i<nr ; i++ ) {
01490       som[i] += 2*xci[i] ;
01491       xco[i] = som[i] ;
01492       som[i] = -som[i] ;
01493     }   // Fin de la boucle sur r
01494       }   // Fin de la boucle sur theta
01495       // Normalisation du premier theta
01496       for (int i=0 ; i<nr ; i++) {
01497     xco[i] *= .5 ;
01498       }
01499       // Positionnement phi suivant
01500       xci += nr*nt ;
01501       xco += nr*nt ;
01502       break ;
01503       
01504     case 0:     // m pair: sin(pair)
01505       // Dernier j: j = nt-1
01506       // Positionnement
01507       xci += nr * nt ;
01508       xco += nr * (nt-1) ;
01509       
01510       // Initialisation de som
01511       for (int i=0 ; i<nr ; i++) {
01512     som[i] = 0. ;
01513     xco[i] = 0. ;
01514       }
01515       
01516       // j suivants
01517       for (int j=nt-2 ; j >= 0 ; j--) {
01518     // Positionnement
01519     xci -= nr ;
01520     xco -= nr ;
01521     // Calcul
01522     for (int i=0 ; i<nr ; i++ ) {
01523       som[i] += 2*xci[i] ;
01524       xco[i] = som[i] ;
01525       som[i] = -som[i] ;
01526     }   // Fin de la boucle sur r
01527       }   // Fin de la boucle sur theta
01528       // Positionnement phi suivant
01529       xci -= nr ;
01530       xci += nr*nt ;
01531       xco += nr*nt ;
01532       break ;
01533     }   // Fin du switch
01534   }   // Fin de la boucle en phi
01535   
01536   // On remet les choses la ou il faut
01537   delete [] tb->t ;
01538   tb->t = xo ;
01539   
01540   //Menage
01541   delete [] som ;
01542   
01543   // base de developpement
01544   int base_r = b & MSQ_R ;
01545   int base_p = b & MSQ_P ;
01546   b = base_r | base_p | T_COSSIN_SI ;
01547     
01548 }
01549 
01550             //----------------------
01551             // cas T_COSSIN_C
01552             //----------------------
01553 void _scost_t_cossin_c(Tbl* tb, int & b) {
01554 
01555     // Peut-etre rien a faire ?
01556     if (tb->get_etat() == ETATZERO) {
01557     int base_r = b & MSQ_R ;
01558     int base_p = b & MSQ_P ;
01559     switch(base_r){
01560         case(R_CHEBPI_P):
01561         b = R_CHEBPI_I | base_p | T_COSSIN_C ;
01562         break ;
01563         case(R_CHEBPI_I):
01564         b = R_CHEBPI_P | base_p | T_COSSIN_C ;
01565         break ;  
01566         default:
01567         b = base_r | base_p | T_COSSIN_C ;
01568         break;
01569     }
01570     return ;
01571     }
01572     
01573     // Protection
01574     assert(tb->get_etat() == ETATQCQ) ;
01575     
01576     // Pour le confort : nbre de points reels.
01577     int nr = (tb->dim).dim[0] ;
01578     int nt = (tb->dim).dim[1] ;
01579     int np = (tb->dim).dim[2] ;
01580     np = np - 2 ;
01581     
01582     // zone de travail
01583     double* somP = new double [nr] ;
01584     double* somN = new double [nr] ;
01585     
01586     // pt. sur le tableau de double resultat
01587     double* xo = new double[(tb->dim).taille] ;
01588     
01589     // Initialisation a zero :
01590     for (int i=0; i<(tb->dim).taille; i++) xo[i] = 0 ; 
01591     
01592     // On y va...
01593     double* xi = tb->t ;
01594     double* xci = xi ;  // Pointeurs
01595     double* xco = xo ;  //  courants
01596     
01597     // k = 0
01598     
01599     // Dernier j: j = nt-1
01600     // Positionnement
01601     xci += nr * (nt-1) ;
01602     xco += nr * (nt-1) ;
01603     
01604     // Initialisation de som 
01605     for (int i=0 ; i<nr ; i++) {
01606     somP[i] = 0. ;
01607     somN[i] = 0. ;
01608     xco[i] = 0. ;   // mise a zero du dernier coef en theta
01609     }
01610     
01611     // j suivants
01612     for (int j=nt-2 ; j >= 0 ; j--) {
01613     // Positionnement
01614     xci -= nr ;
01615     xco -= nr ;
01616     // Calcul
01617     for (int i=0 ; i<nr ; i++ ) {
01618         if((j%2) == 1) {
01619         somN[i] = -somN[i] ;
01620         somN[i] += 2*xci[i] ;
01621         xco[i] = somP[i] ;
01622         }
01623         else {
01624         somP[i] = -somP[i] ;
01625         somP[i] += 2*xci[i] ; 
01626         xco[i] = somN[i] ;
01627         }
01628     }   // Fin de la boucle sur r
01629     }   // Fin de la boucle sur theta
01630     // j=0 : le facteur 2...
01631     for (int i=0 ; i<nr ; i++) xco[i] *= .5 ;
01632  
01633     // Positionnement phi suivant   
01634     xci += nr*nt ;
01635     xco += nr*nt ;
01636     
01637     // k=1
01638     xci += nr*nt ;
01639     xco += nr*nt ;
01640     
01641     // k >= 2
01642     for (int k=2 ; k<np+1 ; k++) {
01643     
01644         // Dernier j: j = nt-1
01645     // Positionnement
01646     xco += nr * (nt-1) ;
01647     xci += nr * (nt-1) ;
01648     
01649     // Initialisation de som
01650     for (int i=0 ; i<nr ; i++) {
01651         somP[i] = 0. ;
01652         somN[i] = 0. ;
01653         xco[i] = 0. ;
01654     }
01655     
01656     // j suivants
01657     for (int j=nt-2 ; j >= 0 ; j--) {
01658         // Positionnement
01659         xci -= nr ;
01660         xco -= nr ;
01661         // Calcul      
01662         for (int i=0 ; i<nr ; i++ ) {
01663         if((j%2) == 1) {
01664             somN[i] = -somN[i];
01665             somN[i] += 2 * xci[i] ;
01666             xco[i] = somP[i] ;
01667         }
01668         else {
01669             somP[i] = -somP[i];
01670             somP[i] += 2 * xci[i] ;
01671             xco[i] = somN[i];
01672         }
01673         }   // Fin de la boucle sur r
01674     }   // Fin de la boucle sur theta
01675     double fac_m = ( (k/2)%2 == 1 ? 0. : 0.5) ;
01676     // j=0 : sin(0*theta) ou facteur 2, suivant la parite de m
01677     for (int i=0 ; i<nr ; i++) xco[i] *= fac_m ;
01678     
01679     // Positionnement phi suivant
01680     xci += nr*nt ;
01681     xco += nr*nt ;
01682     }   // Fin de la boucle sur phi
01683     
01684     // On remet les choses la ou il faut
01685     delete [] tb->t ;
01686     tb->t = xo ;
01687     
01688     //Menage
01689     delete [] somP ;
01690     delete [] somN ;
01691     
01692     // base de developpement
01693     int base_r = b & MSQ_R ;
01694     int base_p = b & MSQ_P ;
01695     switch(base_r){
01696     case(R_CHEBPI_P):
01697         b = R_CHEBPI_I | base_p | T_COSSIN_C ;
01698         break ;
01699     case(R_CHEBPI_I):
01700         b = R_CHEBPI_P | base_p | T_COSSIN_C ;
01701         break ;  
01702     default:
01703         b = base_r | base_p | T_COSSIN_C ;
01704         break;
01705     }
01706 }
01707 
01708             //--------------------- 
01709             // cas T_COSSIN_S
01710             //----------------
01711 void _scost_t_cossin_s(Tbl* tb, int & b) {
01712 
01713     // Peut-etre rien a faire ?
01714     if (tb->get_etat() == ETATZERO) {
01715     int base_r = b & MSQ_R ;
01716     int base_p = b & MSQ_P ;
01717     switch(base_r){
01718         case(R_CHEBPI_P):
01719         b = R_CHEBPI_I | base_p | T_COSSIN_S ;
01720         break ;
01721         case(R_CHEBPI_I):
01722         b = R_CHEBPI_P | base_p | T_COSSIN_S ;
01723         break ;  
01724         default:
01725         b = base_r | base_p | T_COSSIN_S ;
01726         break;
01727     }
01728     return ;
01729     }
01730     
01731     // Protection
01732     assert(tb->get_etat() == ETATQCQ) ;
01733     
01734     // Pour le confort : nbre de points reels.
01735     int nr = (tb->dim).dim[0] ;
01736     int nt = (tb->dim).dim[1] ;
01737     int np = (tb->dim).dim[2] ;
01738     np = np - 2 ;
01739     
01740     // zone de travail
01741     double* somP = new double [nr] ;
01742     double* somN = new double [nr] ;
01743 
01744     // pt. sur le tableau de double resultat
01745     double* xo = new double[(tb->dim).taille] ;
01746     
01747     // Initialisation a zero :
01748     for (int i=0; i<(tb->dim).taille; i++) xo[i] = 0 ; 
01749     
01750     // On y va...
01751     double* xi = tb->t ;
01752     double* xci = xi ;  // Pointeurs
01753     double* xco = xo ;  //  courants
01754     
01755     // k = 0
01756     
01757     // Dernier j: j = nt-1
01758     // Positionnement
01759     xci += nr * (nt-1) ;
01760     xco += nr * (nt-1) ;
01761     
01762     // Initialisation de som 
01763     for (int i=0 ; i<nr ; i++) {
01764     somP[i] = 0. ;
01765     somN[i] = 0. ;
01766     xco[i] = 0. ;   // mise a zero du dernier coef en theta
01767     }
01768     
01769     // j suivants
01770     for (int j=nt-2 ; j >= 0 ; j--) {
01771     // Positionnement
01772     xci -= nr ;
01773     xco -= nr ;
01774     // Calcul
01775     for (int i=0 ; i<nr ; i++ ) {
01776         if((j%2) == 1) {
01777         somN[i] = -somN[i] ;
01778         somN[i] += 2*xci[i] ;
01779         xco[i] = somP[i] ;
01780         }
01781         else {
01782         somP[i] = -somP[i] ;
01783         somP[i] += 2*xci[i] ;         
01784         xco[i] = somN[i] ;
01785         }
01786     }   // Fin de la boucle sur r
01787     }   // Fin de la boucle sur theta
01788     // j=0 : sin(0*theta)
01789     for (int i=0 ; i<nr ; i++) xco[i] = 0. ;
01790    
01791     // Positionnement phi suivant
01792     xci += nr*nt ;
01793     xco += nr*nt ;
01794     
01795     // k=1
01796     xci += nr*nt ;
01797     xco += nr*nt ;
01798     
01799     // k >= 2
01800     for (int k=2 ; k<np+1 ; k++) {
01801     
01802         // Dernier j: j = nt-1
01803     // Positionnement
01804     xco += nr * (nt-1) ;
01805     xci += nr * (nt-1) ;
01806 
01807     // Initialisation de som
01808     for (int i=0 ; i<nr ; i++) {
01809         somP[i] = 0. ;
01810         somN[i] = 0. ;
01811         xco[i] = 0. ;
01812     }
01813     
01814     // j suivants
01815     for (int j=nt-2 ; j >= 0 ; j--) {
01816         // Positionnement
01817         xci -= nr ;
01818         xco -= nr ;
01819         // Calcul      
01820         for (int i=0 ; i<nr ; i++ ) {
01821         if((j%2) == 1) {
01822           somN[i] = -somN[i] ;
01823           somN[i] += 2*xci[i] ;
01824           xco[i] = somP[i] ;
01825           }
01826           else {
01827           somP[i] = -somP[i] ;
01828           somP[i] += 2*xci[i] ;       
01829           xco[i] = somN[i] ;
01830           }
01831         }   // Fin de la boucle sur r
01832     }   // Fin de la boucle sur theta
01833     
01834     double fac_m = ( (k/2)%2 == 0 ? 0. : 0.5) ;
01835     // j=0 : sin(0*theta) ou facteur 2, suivant la parite de m
01836     for (int i=0 ; i<nr ; i++) xco[i] *= fac_m ;
01837     
01838     // Positionnement phi suivant
01839     xci += nr*nt ;
01840     xco += nr*nt ;
01841   
01842     }   // Fin de la boucle sur phi
01843     
01844     // On remet les choses la ou il faut
01845     delete [] tb->t ;
01846     tb->t = xo ;
01847     
01848     //Menage
01849     delete [] somP ;
01850     delete [] somN ;
01851     
01852     // base de developpement
01853     int base_r = b & MSQ_R ;
01854     int base_p = b & MSQ_P ;
01855     switch(base_r){
01856     case(R_CHEBPI_P):
01857         b = R_CHEBPI_I | base_p | T_COSSIN_S ;
01858         break ;
01859     case(R_CHEBPI_I):
01860         b = R_CHEBPI_P | base_p | T_COSSIN_S ;
01861         break ;  
01862     default:
01863         b = base_r | base_p | T_COSSIN_S ;
01864         break;
01865     }
01866 }
Generated on Tue Feb 7 01:35:18 2012 for LORENE by  doxygen 1.4.6