23char cmp_raccord_C[] =
"$Header: /cvsroot/Lorene/C++/Source/Cmp/cmp_raccord.C,v 1.4 2014/10/13 08:52:48 j_novak Exp $" ;
62Matrice matrice_raccord_pair (
int cont,
double alpha_kernel) {
64 Matrice systeme (cont, cont) ;
65 systeme.set_etat_qcq() ;
66 for (
int i=0 ; i<cont ; i++)
67 for (
int j=0 ; j<cont ; j++)
68 systeme.set(i, j) = 0 ;
71 for (
int i=0 ; i<cont ; i++)
72 for (
int k=0 ; k<cont ; k++)
75 for (
int boucle=0 ; boucle<k ; boucle++)
76 somme *= (4*i*i-boucle*boucle)/(2.*boucle+1.)/alpha_kernel ;
77 systeme.set(k, i) = somme ;
79 int inf = (cont%2 == 1) ? (cont-1)/2 : (cont-2)/2 ;
80 systeme.set_band (cont-1, inf) ;
85Matrice matrice_raccord_impair (
int cont,
double alpha_kernel) {
87 Matrice systeme (cont, cont) ;
88 systeme.set_etat_qcq() ;
89 for (
int i=0 ; i<cont ; i++)
90 for (
int j=0 ; j<cont ; j++)
91 systeme.set(i, j) = 0 ;
94 for (
int i=0 ; i<cont ; i++)
95 for (
int k=0 ; k<cont ; k++)
98 for (
int boucle=0 ; boucle<k ; boucle++)
99 somme *= (
pow(2*i+1, 2.)-boucle*boucle)/(2.*boucle+1.)/alpha_kernel ;
100 systeme.set(k, i) = somme ;
102 int inf = (cont%2 == 0) ? cont/2 : (cont-1)/2 ;
103 systeme.set_band (cont-1, inf) ;
109Tbl sec_membre_raccord (Tbl coef,
int cont,
double alpha_shell) {
111 assert (coef.get_etat() != ETATNONDEF) ;
112 int nr = coef.get_dim(0) ;
114 Tbl sec_membre(cont) ;
115 sec_membre.set_etat_qcq() ;
116 for (
int i=0 ; i<cont ; i++)
117 sec_membre.set(i) = 0 ;
120 for (
int i=0 ; i<nr ; i++)
121 for (
int k=0 ; k<cont ; k++)
124 for (
int boucle=0 ; boucle<k ; boucle++)
125 somme *= (i*i-boucle*boucle)/(2.*boucle+1.)/alpha_shell ;
127 sec_membre.set(k) += coef(i)*somme ;
129 sec_membre.set(k) -= coef(i)*somme ;
136Tbl regularise (Tbl coef,
int nr,
int base_r) {
139 assert (coef.get_etat() != ETATNONDEF) ;
140 int cont = coef.get_dim(0) ;
141 assert (nr >= cont) ;
144 resu.set_etat_qcq() ;
146 double* x4coef =
new double[nr] ;
147 for (
int i=0 ; i<cont ; i++)
148 x4coef[i] = coef(i) ;
149 for (
int i=cont ; i<nr ; i++)
151 double* x6coef =
new double[nr] ;
153 multx2_1d (nr, &x4coef, base_r) ;
154 multx2_1d (nr, &x4coef, base_r) ;
155 for (
int i=0 ; i<nr ; i++)
156 x6coef[i] = x4coef[i] ;
157 multx2_1d (nr, &x6coef, base_r) ;
159 for (
int i=0 ; i<nr ; i++)
160 resu.set(i) = 3*x4coef[i]-2*x6coef[i] ;
171 assert (
etat != ETATNONDEF) ;
178 if (mapping == 0x0) {
180 "raccord : The mapping does not belong to the class Map_af !"
205 double alpha_kernel = mapping->
get_alpha()[0] ;
206 double alpha_shell = mapping->
get_alpha()[1] ;
208 int base_r, m_quant, l_quant ;
210 for (
int k=0 ; k<nbrep_kernel+1 ; k++)
211 for (
int j=0 ; j<nbret_kernel ; j++)
212 if (nullite_plm(j, nbret_kernel, k,nbrep_kernel,
va.
base) == 1)
213 if (nullite_plm(j, nbret_shell, k, nbrep_shell,
va.
base) == 1)
216 donne_lm(nz, 0, j, k,
va.
base, m_quant, l_quant, base_r) ;
221 Tbl facteur (nbrer_kernel) ;
223 for (
int i=0 ; i<nbrer_shell ; i++)
225 facteur.
set(i) = (*
va.
c_cf)(1, k, j, i) ;
227 Tbl sec_membre (sec_membre_raccord (facteur, cont, alpha_shell)) ;
230 systeme = matrice_raccord_pair (cont, alpha_kernel) ;
232 systeme = matrice_raccord_impair (cont, alpha_kernel) ;
235 Tbl regulier (nbrer_kernel) ;
238 for (
int i=0 ; i<cont ; i++)
242 regulier = regularise (soluce, nbrer_kernel,
R_CHEBP) ;
244 regulier = regularise (soluce, nbrer_kernel,
R_CHEBI) ;
246 for (
int i=0 ; i<nbrer_kernel ; i++)
Valeur va
The numerical value of the Cmp
int etat
Logical state (ETATNONDEF , ETATQCQ or ETATZERO ).
const Map * get_mp() const
Returns the mapping.
void raccord(int n)
Performs the matching of the nucleus with respect to the first shell.
const double * get_alpha() const
Returns the pointer on the array alpha.
const Mg3d * get_mg() const
Gives the Mg3d on which the mapping is defined.
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_nt(int l) const
Returns the number of points in the co-latitude direction ( ) in domain no. l.
int get_nzone() const
Returns the number of domains.
int get_nr(int l) const
Returns the number of points in the radial direction ( ) in domain no. l.
int get_type_r(int l) const
Returns the type of sampling in the radial direction in domain no.
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
Tbl ** t
Array (size nzone ) of pointers on the Tbl 's which contain the spectral coefficients in each domain.
void annule_hard()
Sets the Tbl to zero in a hard way.
double & set(int i)
Read/write of a particular element (index i) (1D case)
void set_etat_cf_qcq()
Sets the logical state to ETATQCQ (ordinary state) for values in the configuration space (Mtbl_cf c_c...
void ylm()
Computes the coefficients of *this.
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void coef() const
Computes the coeffcients of *this.
void ylm_i()
Inverse of ylm()
Base_val base
Bases on which the spectral expansion is performed.
Cmp pow(const Cmp &, int)
Power .
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEBP
base de Cheb. paire (rare) seulement