op_ssint.C

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