64 #include "utilitaires.h"
67 #include "param_elliptic.h"
70 #include "sym_tensor.h"
86 assert(mp_aff != 0x0) ;
108 assert(aaa.
get_etat() != ETATNONDEF) ;
111 bool ced = (mgrid.
get_type_r(nzm1) == UNSURR) ;
112 int n_shell = ced ? nz-2 : nzm1 ;
116 n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
117 bool cedbc = (ced && (nz_bc == nzm1)) ;
120 assert(par_bc != 0x0) ;
124 int nt = mgrid.
get_nt(0) ;
125 int np = mgrid.
get_np(0) ;
153 int l_q, m_q, base_r ;
162 int nr = mgrid.
get_nr(lz) ;
167 for (
int k=0 ; k<np+1 ; k++) {
168 for (
int j=0 ; j<nt ; j++) {
171 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
176 for (
int i=0; i<nr; i++) {
177 sol_part_mu.
set(lz, k, j, i) = 0 ;
178 sol_part_x.
set(lz, k, j, i) = 0 ;
180 for (
int i=0; i<nr; i++) {
181 sol_hom2_mu.
set(lz, k, j, i) = 0 ;
182 sol_hom2_x.
set(lz, k, j, i) = 0 ;
192 for (
int lz=1; lz <= n_shell; lz++) {
193 int nr = mgrid.
get_nr(lz) ;
194 assert(mgrid.
get_nt(lz) == nt) ;
195 assert(mgrid.
get_np(lz) == np) ;
197 double ech = mp_aff->
get_beta()[lz] / alpha ;
201 for (
int k=0 ; k<np+1 ; k++) {
202 for (
int j=0 ; j<nt ; j++) {
205 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
210 for (
int lin=0; lin<nr; lin++)
211 for (
int col=0; col<nr; col++)
212 ope.
set(lin,col) = mxd(lin,col) + ech*md(lin,col)
214 for (
int lin=0; lin<nr; lin++)
215 for (
int col=0; col<nr; col++)
216 ope.
set(lin,col+nr) = (2-l_q*(l_q+1))*mid(lin,col) ;
217 for (
int lin=0; lin<nr; lin++)
218 for (
int col=0; col<nr; col++)
219 ope.
set(lin+nr,col) = -mid(lin,col) ;
220 for (
int lin=0; lin<nr; lin++)
221 for (
int col=0; col<nr; col++)
222 ope.
set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col) ;
226 for (
int col=0; col<2*nr; col++) {
227 ope.
set(ind0+nr-1, col) = 0 ;
228 ope.
set(ind1+nr-1, col) = 0 ;
230 ope.
set(ind0+nr-1, ind0) = 1 ;
231 ope.
set(ind1+nr-1, ind1) = 1 ;
237 for (
int lin=0; lin<nr; lin++)
239 for (
int lin=0; lin<nr; lin++)
242 sec.
set(ind0+nr-1) = 0 ;
243 sec.
set(ind1+nr-1) = 0 ;
245 for (
int i=0; i<nr; i++) {
246 sol_part_mu.
set(lz, k, j, i) = sol(i) ;
247 sol_part_x.
set(lz, k, j, i) = sol(i+nr) ;
250 sec.
set(ind0+nr-1) = 1 ;
252 for (
int i=0; i<nr; i++) {
253 sol_hom1_mu.
set(lz, k, j, i) = sol(i) ;
254 sol_hom1_x.
set(lz, k, j, i) = sol(i+nr) ;
256 sec.
set(ind0+nr-1) = 0 ;
257 sec.
set(ind1+nr-1) = 1 ;
259 for (
int i=0; i<nr; i++) {
260 sol_hom2_mu.
set(lz, k, j, i) = sol(i) ;
261 sol_hom2_x.
set(lz, k, j, i) = sol(i+nr) ;
271 if (cedbc) {
int lz = nzm1 ;
272 int nr = mgrid.
get_nr(lz) ;
273 assert(mgrid.
get_nt(lz) == nt) ;
274 assert(mgrid.
get_np(lz) == np) ;
279 for (
int k=0 ; k<np+1 ; k++) {
280 for (
int j=0 ; j<nt ; j++) {
283 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
287 for (
int lin=0; lin<nr; lin++)
288 for (
int col=0; col<nr; col++)
289 ope.
set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
290 for (
int lin=0; lin<nr; lin++)
291 for (
int col=0; col<nr; col++)
292 ope.
set(lin,col+nr) = (2-l_q*(l_q+1))*ms(lin,col) ;
293 for (
int lin=0; lin<nr; lin++)
294 for (
int col=0; col<nr; col++)
295 ope.
set(lin+nr,col) = -ms(lin,col) ;
296 for (
int lin=0; lin<nr; lin++)
297 for (
int col=0; col<nr; col++)
298 ope.
set(lin+nr,col+nr) = -md(lin,col) ;
303 for (
int col=0; col<2*nr; col++) {
304 ope.
set(ind0+nr-1, col) = 0 ;
305 ope.
set(ind1+nr-2, col) = 0 ;
306 ope.
set(ind1+nr-1, col) = 0 ;
308 for (
int col=0; col<nr; col++) {
309 ope.
set(ind0+nr-1, col+ind0) = 1 ;
310 ope.
set(ind1+nr-1, col+ind1) = 1 ;
312 ope.
set(ind1+nr-2, ind1+1) = 1 ;
318 for (
int lin=0; lin<nr; lin++)
320 for (
int lin=0; lin<nr; lin++)
323 sec.
set(ind0+nr-1) = 0 ;
324 sec.
set(ind1+nr-2) = 0 ;
325 sec.
set(ind1+nr-1) = 0 ;
327 for (
int i=0; i<nr; i++) {
328 sol_part_mu.
set(lz, k, j, i) = sol(i) ;
329 sol_part_x.
set(lz, k, j, i) = sol(i+nr) ;
332 sec.
set(ind1+nr-2) = 1 ;
334 for (
int i=0; i<nr; i++) {
335 sol_hom1_mu.
set(lz, k, j, i) = sol(i) ;
336 sol_hom1_x.
set(lz, k, j, i) = sol(i+nr) ;
346 int taille = 2*nz_bc ;
347 if (cedbc) taille-- ;
351 Tbl sec_membre(taille) ;
352 Matrice systeme(taille, taille) ;
353 int ligne ;
int colonne ;
356 double c_mu = (cedbc ? 0 : par_bc->
get_tbl_mod(0)(0) ) ;
357 double d_mu = (cedbc ? 0 : par_bc->
get_tbl_mod(0)(1) ) ;
358 double c_x = (cedbc ? 0 : par_bc->
get_tbl_mod(0)(2) ) ;
359 double d_x = (cedbc ? 0 : par_bc->
get_tbl_mod(0)(3) ) ;
360 Mtbl_cf dhom1_mu = sol_hom1_mu ;
361 Mtbl_cf dhom2_mu = sol_hom2_mu ;
362 Mtbl_cf dpart_mu = sol_part_mu ;
378 for (
int k=0 ; k<np+1 ; k++)
379 for (
int j=0 ; j<nt ; j++) {
381 if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
389 int nr = mgrid.
get_nr(1) ;
427 systeme.
set(ligne, colonne) =
429 systeme.
set(ligne, colonne+1) =
435 systeme.
set(ligne, colonne) =
437 systeme.
set(ligne, colonne+1) =
445 for (
int zone=2 ; zone<nz_bc ; zone++) {
450 systeme.
set(ligne, colonne) =
452 systeme.
set(ligne, colonne+1) =
458 systeme.
set(ligne, colonne) =
460 systeme.
set(ligne, colonne+1) =
467 systeme.
set(ligne, colonne) =
469 systeme.
set(ligne, colonne+1) =
475 systeme.
set(ligne, colonne) =
477 systeme.
set(ligne, colonne+1) =
486 nr = mgrid.
get_nr(nz_bc) ;
487 double alpha = mp_aff->
get_alpha()[nz_bc] ;
491 systeme.
set(ligne, colonne) =
493 if (!cedbc) systeme.
set(ligne, colonne+1) =
499 systeme.
set(ligne, colonne) =
501 if (!cedbc) systeme.
set(ligne, colonne+1) =
509 systeme.
set(ligne, colonne) =
514 systeme.
set(ligne, colonne+1) =
543 for (
int i=0 ; i<nr ; i++) {
544 mmu.
set(0, k, j, i) = 0 ;
545 mw.
set(0, k, j, i) = 0 ;
548 for (
int i=0 ; i<nr ; i++) {
549 mmu.
set(1, k, j, i) = sol_part_mu(1, k, j, i)
550 + facteur(conte)*sol_hom1_mu(1, k, j, i)
551 + facteur(conte+1)*sol_hom2_mu(1, k, j, i);
552 mw.
set(1, k, j, i) = sol_part_x(1, k, j, i)
553 + facteur(conte)*sol_hom1_x(1, k, j, i)
554 + facteur(conte+1)*sol_hom2_x(1,k,j,i);
557 for (
int zone=2 ; zone<=n_shell ; zone++) {
559 for (
int i=0 ; i<nr ; i++) {
560 mmu.
set(zone, k, j, i) = sol_part_mu(zone, k, j, i)
561 + facteur(conte)*sol_hom1_mu(zone, k, j, i)
562 + facteur(conte+1)*sol_hom2_mu(zone, k, j, i) ;
564 mw.
set(zone, k, j, i) = sol_part_x(zone, k, j, i)
565 + facteur(conte)*sol_hom1_x(zone, k, j, i)
566 + facteur(conte+1)*sol_hom2_x(zone, k, j, i) ;
572 for (
int i=0 ; i<nr ; i++) {
573 mmu.
set(nzm1, k, j, i) = sol_part_mu(nzm1, k, j, i)
574 + facteur(conte)*sol_hom1_mu(nzm1, k, j, i) ;
576 mw.
set(nzm1, k, j, i) = sol_part_x(nzm1, k, j, i)
577 + facteur(conte)*sol_hom1_x(nzm1, k, j, i) ;
607 assert(mp_aff != 0x0) ;
614 bool ced = (mgrid.
get_type_r(nzm1) == UNSURR) ;
615 int n_shell = ced ? nz-2 : nzm1 ;
620 n_shell = (nz_bc < n_shell ? nz_bc : n_shell) ;
621 bool cedbc = (ced && (nz_bc == nzm1)) ;
624 assert(par_bc != 0x0) ;
628 int nt = mgrid.
get_nt(0) ;
629 int np = mgrid.
get_np(0) ;
631 int l_q, m_q, base_r;
681 assert (bbt.
get_etat() != ETATNONDEF) ;
682 assert (hh.
get_etat() != ETATNONDEF) ;
691 bool bnull = (bbt.
get_etat() == ETATZERO) ;
704 bool hnull = (hh.
get_etat() == ETATZERO) ;
711 bool need_calculation = true ;
712 if (par_mat != 0x0) {
713 bool param_new = false ;
718 if (par_mat->
get_int_mod(0) < nz_bc) param_new =
true ;
719 if (par_mat->
get_int_mod(1) != lmax) param_new =
true ;
725 for (
int l=1; l<= n_shell; l++) {
728 mp_aff->
get_alpha()[l]) > 2.e-15) param_new = true ;
741 int* nz_bc_new =
new int(nz_bc) ;
743 int* lmax_new =
new int(lmax) ;
745 int* type_t_new =
new int(mgrid.
get_type_t()) ;
747 int* type_p_new =
new int(mgrid.
get_type_p()) ;
752 for (
int l=0; l<nz; l++)
754 Tbl* palpha =
new Tbl(nz) ;
758 for (
int l=1; l<nzm1; l++)
762 else need_calculation = false ;
778 sol_Dirac_l01_bound(hh2, hrr, tilde_eta, bound_hrr, bound_eta, par_mat) ;
794 Itbl mat_done(lmax) ;
804 int nr = mgrid.
get_nr(lz) ;
808 if (need_calculation && (par_mat != 0x0)) mat_done.
annule_hard() ;
810 for (
int k=0 ; k<np+1 ; k++) {
811 for (
int j=0 ; j<nt ; j++) {
814 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
815 if (need_calculation) {
820 for (
int lin=0; lin<nr; lin++)
821 for (
int col=0; col<nr; col++)
822 ope.
set(lin,col) = md(lin,col) + 3*ms(lin,col) ;
823 for (
int lin=0; lin<nr; lin++)
824 for (
int col=0; col<nr; col++)
825 ope.
set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
826 for (
int lin=0; lin<nr; lin++)
827 for (
int col=0; col<nr; col++)
828 ope.
set(lin,col+2*nr) = 0 ;
829 for (
int lin=0; lin<nr; lin++)
830 for (
int col=0; col<nr; col++)
831 ope.
set(lin+nr,col) = -0.5*ms(lin,col) ;
832 for (
int lin=0; lin<nr; lin++)
833 for (
int col=0; col<nr; col++)
834 ope.
set(lin+nr,col+nr) = md(lin,col) + 3*ms(lin,col) ;
835 for (
int lin=0; lin<nr; lin++)
836 for (
int col=0; col<nr; col++)
837 ope.
set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*ms(lin,col) ;
838 for (
int lin=0; lin<nr; lin++)
839 for (
int col=0; col<nr; col++)
840 ope.
set(lin+2*nr,col) = -0.5*md(lin,col)/double(l_q+1)
841 - 0.5*double(l_q+4)/double(l_q+1)*ms(lin,col) ;
842 for (
int lin=0; lin<nr; lin++)
843 for (
int col=0; col<nr; col++)
844 ope.
set(lin+2*nr,col+nr) = -2*ms(lin,col) ;
845 for (
int lin=0; lin<nr; lin++)
846 for (
int col=0; col<nr; col++)
847 ope.
set(lin+2*nr,col+2*nr) = (l_q+2)*md(lin,col)
848 + l_q*(l_q+2)*ms(lin,col) ;
851 for (
int col=0; col<3*nr; col++)
852 if (l_q>2) ope.
set(ind2+nr-2, col) = 0 ;
853 for (
int col=0; col<3*nr; col++) {
854 ope.
set(nr-1, col) = 0 ;
855 ope.
set(2*nr-1, col) = 0 ;
856 ope.
set(3*nr-1, col) = 0 ;
860 for (
int col=0; col<nr; col++) {
861 ope.
set(nr-1, col) = pari ;
862 ope.
set(2*nr-1, col+nr) = pari ;
863 ope.
set(3*nr-1, col+2*nr) = pari ;
868 ope.
set(nr-1, nr-1) = 1 ;
869 ope.
set(2*nr-1, 2*nr-1) = 1 ;
870 ope.
set(3*nr-1, 3*nr-1) = 1 ;
873 ope.
set(ind2+nr-2, ind2) = 1 ;
876 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
879 mat_done.
set(l_q) = 1 ;
883 const Matrice& oper = (par_mat == 0x0 ? ope :
888 for (
int lin=0; lin<2*nr; lin++)
890 for (
int lin=0; lin<nr; lin++)
895 for (
int lin=0; lin<nr; lin++)
897 for (
int lin=0; lin<nr; lin++)
901 for (
int lin=0; lin<nr; lin++)
902 sec.
set(2*nr+lin)= -0.5/double(l_q+1)*(
907 for (
int lin=0; lin<nr; lin++)
908 sec.
set(2*nr+lin) = -0.5/double(l_q+1)*(
914 if (l_q>2) sec.
set(ind2+nr-2) = 0 ;
915 sec.
set(3*nr-1) = 0 ;
917 for (
int i=0; i<nr; i++) {
918 sol_part_hrr.
set(lz, k, j, i) = sol(i) ;
919 sol_part_eta.
set(lz, k, j, i) = sol(i+nr) ;
920 sol_part_w.
set(lz, k, j, i) = sol(i+2*nr) ;
924 sec.
set(ind2+nr-2) = 1 ;
933 for (
int i=0; i<nr; i++) {
934 sol_hom3_hrr.
set(lz, k, j, i) = sol(i) ;
935 sol_hom3_eta.
set(lz, k, j, i) = sol(i+nr) ;
936 sol_hom3_w.
set(lz, k, j, i) = sol(i+2*nr) ;
948 for (
int lz=1; lz<= n_shell; lz++) {
949 if (need_calculation && (par_mat != 0x0)) mat_done.
annule_hard() ;
950 int nr = mgrid.
get_nr(lz) ;
955 double ech = mp_aff->
get_beta()[lz] / alpha ;
958 for (
int k=0 ; k<np+1 ; k++) {
959 for (
int j=0 ; j<nt ; j++) {
962 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
963 if (need_calculation) {
969 for (
int lin=0; lin<nr; lin++)
970 for (
int col=0; col<nr; col++)
971 ope.
set(lin,col) = mxd(lin,col) + ech*md(lin,col)
973 for (
int lin=0; lin<nr; lin++)
974 for (
int col=0; col<nr; col++)
975 ope.
set(lin,col+nr) = -l_q*(l_q+1)*mid(lin,col) ;
976 for (
int lin=0; lin<nr; lin++)
977 for (
int col=0; col<nr; col++)
978 ope.
set(lin,col+2*nr) = 0 ;
979 for (
int lin=0; lin<nr; lin++)
980 for (
int col=0; col<nr; col++)
981 ope.
set(lin+nr,col) = -0.5*mid(lin,col) ;
982 for (
int lin=0; lin<nr; lin++)
983 for (
int col=0; col<nr; col++)
984 ope.
set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col)
986 for (
int lin=0; lin<nr; lin++)
987 for (
int col=0; col<nr; col++)
988 ope.
set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*mid(lin,col) ;
989 for (
int lin=0; lin<nr; lin++)
990 for (
int col=0; col<nr; col++)
991 ope.
set(lin+2*nr,col) =
992 -0.5/double(l_q+1)*(mxd(lin,col) + ech*md(lin,col)
993 + double(l_q+4)*mid(lin,col)) ;
994 for (
int lin=0; lin<nr; lin++)
995 for (
int col=0; col<nr; col++)
996 ope.
set(lin+2*nr,col+nr) = -2*mid(lin,col) ;
997 for (
int lin=0; lin<nr; lin++)
998 for (
int col=0; col<nr; col++)
999 ope.
set(lin+2*nr,col+2*nr) =
1000 double(l_q+2)*(mxd(lin,col) + ech*md(lin,col)
1001 + l_q*mid(lin,col)) ;
1002 for (
int col=0; col<3*nr; col++) {
1003 ope.
set(ind0+nr-1, col) = 0 ;
1004 ope.
set(ind1+nr-1, col) = 0 ;
1005 ope.
set(ind2+nr-1, col) = 0 ;
1007 ope.
set(ind0+nr-1, ind0) = 1 ;
1008 ope.
set(ind1+nr-1, ind1) = 1 ;
1009 ope.
set(ind2+nr-1, ind2) = 1 ;
1012 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1015 mat_done.
set(l_q) = 1 ;
1018 const Matrice& oper = (par_mat == 0x0 ? ope :
1023 for (
int lin=0; lin<2*nr; lin++)
1025 for (
int lin=0; lin<nr; lin++)
1030 for (
int lin=0; lin<nr; lin++)
1032 for (
int lin=0; lin<nr; lin++)
1036 for (
int lin=0; lin<nr; lin++)
1037 sec.
set(2*nr+lin) = -0.5/double(l_q+1)*(
1042 for (
int lin=0; lin<nr; lin++)
1043 sec.
set(2*nr+lin) = -0.5/double(l_q+1)*(
1049 sec.
set(ind0+nr-1) = 0 ;
1050 sec.
set(ind1+nr-1) = 0 ;
1051 sec.
set(ind2+nr-1) = 0 ;
1053 for (
int i=0; i<nr; i++) {
1054 sol_part_hrr.
set(lz, k, j, i) = sol(i) ;
1055 sol_part_eta.
set(lz, k, j, i) = sol(i+nr) ;
1056 sol_part_w.
set(lz, k, j, i) = sol(i+2*nr) ;
1059 sec.
set(ind0+nr-1) = 1 ;
1061 for (
int i=0; i<nr; i++) {
1062 sol_hom1_hrr.
set(lz, k, j, i) = sol(i) ;
1063 sol_hom1_eta.
set(lz, k, j, i) = sol(i+nr) ;
1064 sol_hom1_w.
set(lz, k, j, i) = sol(i+2*nr) ;
1066 sec.
set(ind0+nr-1) = 0 ;
1067 sec.
set(ind1+nr-1) = 1 ;
1069 for (
int i=0; i<nr; i++) {
1070 sol_hom2_hrr.
set(lz, k, j, i) = sol(i) ;
1071 sol_hom2_eta.
set(lz, k, j, i) = sol(i+nr) ;
1072 sol_hom2_w.
set(lz, k, j, i) = sol(i+2*nr) ;
1074 sec.
set(ind1+nr-1) = 0 ;
1075 sec.
set(ind2+nr-1) = 1 ;
1077 for (
int i=0; i<nr; i++) {
1078 sol_hom3_hrr.
set(lz, k, j, i) = sol(i) ;
1079 sol_hom3_eta.
set(lz, k, j, i) = sol(i+nr) ;
1080 sol_hom3_w.
set(lz, k, j, i) = sol(i+2*nr) ;
1090 if (cedbc) {
int lz = nzm1 ;
1091 if (need_calculation && (par_mat != 0x0)) mat_done.
annule_hard() ;
1092 int nr = mgrid.
get_nr(lz) ;
1096 double alpha = mp_aff->
get_alpha()[lz] ;
1099 for (
int k=0 ; k<np+1 ; k++) {
1100 for (
int j=0 ; j<nt ; j++) {
1103 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
1104 if (need_calculation) {
1109 for (
int lin=0; lin<nr; lin++)
1110 for (
int col=0; col<nr; col++)
1111 ope.
set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
1112 for (
int lin=0; lin<nr; lin++)
1113 for (
int col=0; col<nr; col++)
1114 ope.
set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
1115 for (
int lin=0; lin<nr; lin++)
1116 for (
int col=0; col<nr; col++)
1117 ope.
set(lin,col+2*nr) = 0 ;
1118 for (
int lin=0; lin<nr; lin++)
1119 for (
int col=0; col<nr; col++)
1120 ope.
set(lin+nr,col) = -0.5*ms(lin,col) ;
1121 for (
int lin=0; lin<nr; lin++)
1122 for (
int col=0; col<nr; col++)
1123 ope.
set(lin+nr,col+nr) = -md(lin,col) + 3*ms(lin,col) ;
1124 for (
int lin=0; lin<nr; lin++)
1125 for (
int col=0; col<nr; col++)
1126 ope.
set(lin+nr,col+2*nr) = (2. - l_q*(l_q+1))*ms(lin,col) ;
1127 for (
int lin=0; lin<nr; lin++)
1128 for (
int col=0; col<nr; col++)
1129 ope.
set(lin+2*nr,col) = 0.5*md(lin,col)/double(l_q+1)
1130 - 0.5*double(l_q+4)/double(l_q+1)*ms(lin,col) ;
1131 for (
int lin=0; lin<nr; lin++)
1132 for (
int col=0; col<nr; col++)
1133 ope.
set(lin+2*nr,col+nr) = -2*ms(lin,col) ;
1134 for (
int lin=0; lin<nr; lin++)
1135 for (
int col=0; col<nr; col++)
1136 ope.
set(lin+2*nr,col+2*nr) = -(l_q+2)*md(lin,col)
1137 + l_q*(l_q+2)*ms(lin,col) ;
1139 for (
int col=0; col<3*nr; col++) {
1140 ope.
set(ind0+nr-2, col) = 0 ;
1141 ope.
set(ind0+nr-1, col) = 0 ;
1142 ope.
set(ind1+nr-2, col) = 0 ;
1143 ope.
set(ind1+nr-1, col) = 0 ;
1144 ope.
set(ind2+nr-1, col) = 0 ;
1146 for (
int col=0; col<nr; col++) {
1147 ope.
set(ind0+nr-1, col+ind0) = 1 ;
1148 ope.
set(ind1+nr-1, col+ind1) = 1 ;
1149 ope.
set(ind2+nr-1, col+ind2) = 1 ;
1151 ope.
set(ind0+nr-2, ind0+1) = 1 ;
1152 ope.
set(ind1+nr-2, ind1+2) = 1 ;
1155 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1158 mat_done.
set(l_q) = 1 ;
1161 const Matrice& oper = (par_mat == 0x0 ? ope :
1167 for (
int lin=0; lin<2*nr; lin++)
1169 for (
int lin=0; lin<nr; lin++)
1174 for (
int lin=0; lin<nr; lin++)
1176 for (
int lin=0; lin<nr; lin++)
1180 for (
int lin=0; lin<nr; lin++)
1181 sec.
set(2*nr+lin) = -0.5/double(l_q+1)*(
1186 for (
int lin=0; lin<nr; lin++)
1187 sec.
set(2*nr+lin) = -0.5/double(l_q+1)*(
1193 sec.
set(ind0+nr-2) = 0 ;
1194 sec.
set(ind0+nr-1) = 0 ;
1195 sec.
set(ind1+nr-1) = 0 ;
1196 sec.
set(ind1+nr-2) = 0 ;
1197 sec.
set(ind2+nr-1) = 0 ;
1199 for (
int i=0; i<nr; i++) {
1200 sol_part_hrr.
set(lz, k, j, i) = sol(i) ;
1201 sol_part_eta.
set(lz, k, j, i) = sol(i+nr) ;
1202 sol_part_w.
set(lz, k, j, i) = sol(i+2*nr) ;
1205 sec.
set(ind0+nr-2) = 1 ;
1207 for (
int i=0; i<nr; i++) {
1208 sol_hom1_hrr.
set(lz, k, j, i) = sol(i) ;
1209 sol_hom1_eta.
set(lz, k, j, i) = sol(i+nr) ;
1210 sol_hom1_w.
set(lz, k, j, i) = sol(i+2*nr) ;
1212 sec.
set(ind0+nr-2) = 0 ;
1213 sec.
set(ind1+nr-2) = 1 ;
1215 for (
int i=0; i<nr; i++) {
1216 sol_hom2_hrr.
set(lz, k, j, i) = sol(i) ;
1217 sol_hom2_eta.
set(lz, k, j, i) = sol(i+nr) ;
1218 sol_hom2_w.
set(lz, k, j, i) = sol(i+2*nr) ;
1229 int taille = 3*nz_bc ;
1230 if (cedbc) taille-- ;
1234 Tbl sec_membre(taille) ;
1235 Matrice systeme(taille, taille) ;
1236 int ligne ;
int colonne ;
1239 double chrr = (cedbc ? 0 : par_bc->
get_tbl_mod()(4) ) ;
1240 double dhrr = (cedbc ? 0 : par_bc->
get_tbl_mod()(5) ) ;
1241 double ceta = (cedbc ? 0 : par_bc->
get_tbl_mod()(6) ) ;
1242 double deta = (cedbc ? 0 : par_bc->
get_tbl_mod()(7) ) ;
1243 double cw = (cedbc ? 0 : par_bc->
get_tbl_mod()(8) ) ;
1244 double dw = (cedbc ? 0 : par_bc->
get_tbl_mod()(9) ) ;
1245 Mtbl_cf dhom1_hrr = sol_hom1_hrr ;
1246 Mtbl_cf dhom2_hrr = sol_hom2_hrr ;
1247 Mtbl_cf dhom3_hrr = sol_hom3_hrr ;
1248 Mtbl_cf dpart_hrr = sol_part_hrr ;
1249 Mtbl_cf dhom1_eta = sol_hom1_eta ;
1250 Mtbl_cf dhom2_eta = sol_hom2_eta ;
1251 Mtbl_cf dhom3_eta = sol_hom3_eta ;
1252 Mtbl_cf dpart_eta = sol_part_eta ;
1253 Mtbl_cf dhom1_w = sol_hom1_w ;
1254 Mtbl_cf dhom2_w = sol_hom2_w ;
1255 Mtbl_cf dhom3_w = sol_hom3_w ;
1256 Mtbl_cf dpart_w = sol_part_w ;
1264 Mtbl_cf d2hom1_hrr = dhom1_hrr ;
1265 Mtbl_cf d2hom2_hrr = dhom2_hrr ;
1266 Mtbl_cf d2hom3_hrr = dhom3_hrr;
1267 Mtbl_cf d2part_hrr = dpart_hrr ;
1268 d2hom1_hrr.
dsdx(); d2hom2_hrr.
dsdx(); d2hom3_hrr.
dsdx(); d2part_hrr.
dsdx();
1276 for (
int k=0 ; k<np+1 ; k++)
1277 for (
int j=0 ; j<nt ; j++) {
1279 if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q > 1)) {
1287 int nr = mgrid.
get_nr(1);
1288 double alpha2 = mp_aff->
get_alpha()[1] ;
1307 sec_membre.
set(ligne) += blob2;
1328 systeme.
set(ligne, colonne) =
1330 systeme.
set(ligne, colonne+1) =
1332 systeme.
set(ligne, colonne+2) =
1338 systeme.
set(ligne, colonne) =
1340 systeme.
set(ligne, colonne+1) =
1342 systeme.
set(ligne, colonne+2) =
1348 systeme.
set(ligne, colonne) =
1350 systeme.
set(ligne, colonne+1) =
1352 systeme.
set(ligne, colonne+2) =
1361 for (
int zone=2 ; zone<nz_bc ; zone++) {
1362 nr = mgrid.
get_nr(zone) ;
1366 systeme.
set(ligne, colonne) =
1368 systeme.
set(ligne, colonne+1) =
1370 systeme.
set(ligne, colonne+2) =
1376 systeme.
set(ligne, colonne) =
1378 systeme.
set(ligne, colonne+1) =
1380 systeme.
set(ligne, colonne+2) =
1386 systeme.
set(ligne, colonne) =
1388 systeme.
set(ligne, colonne+1) =
1390 systeme.
set(ligne, colonne+2) =
1397 systeme.
set(ligne, colonne) =
1399 systeme.
set(ligne, colonne+1) =
1401 systeme.
set(ligne, colonne+2) =
1407 systeme.
set(ligne, colonne) =
1409 systeme.
set(ligne, colonne+1) =
1411 systeme.
set(ligne, colonne+2) =
1417 systeme.
set(ligne, colonne) =
1419 systeme.
set(ligne, colonne+1) =
1421 systeme.
set(ligne, colonne+2) =
1430 nr = mgrid.
get_nr(nz_bc) ;
1431 double alpha = mp_aff->
get_alpha()[nz_bc] ;
1435 systeme.
set(ligne, colonne) =
1437 systeme.
set(ligne, colonne+1) =
1439 if (!cedbc) systeme.
set(ligne, colonne+2) =
1445 systeme.
set(ligne, colonne) =
1447 systeme.
set(ligne, colonne+1) =
1449 if (!cedbc) systeme.
set(ligne, colonne+2) =
1455 systeme.
set(ligne, colonne) =
1457 systeme.
set(ligne, colonne+1) =
1459 if (!cedbc) systeme.
set(ligne, colonne+2) =
1466 systeme.
set(ligne, colonne) =
1473 systeme.
set(ligne, colonne+1) =
1480 systeme.
set(ligne, colonne+2) =
1514 for (
int i=0 ; i<nr ; i++) {
1515 mhrr.
set(0, k, j, i) = 0 ;
1516 meta.
set(0, k, j, i) = 0 ;
1517 mw.
set(0, k, j, i) = 0 ;
1519 for (
int zone=1 ; zone<=n_shell ; zone++) {
1520 nr = mgrid.
get_nr(zone) ;
1521 for (
int i=0 ; i<nr ; i++) {
1522 mhrr.
set(zone, k, j, i) = sol_part_hrr(zone, k, j, i)
1523 + facteur(conte)*sol_hom1_hrr(zone, k, j, i)
1524 + facteur(conte+1)*sol_hom2_hrr(zone, k, j, i)
1525 + facteur(conte+2)*sol_hom3_hrr(zone, k, j, i) ;
1527 meta.
set(zone, k, j, i) = sol_part_eta(zone, k, j, i)
1528 + facteur(conte)*sol_hom1_eta(zone, k, j, i)
1529 + facteur(conte+1)*sol_hom2_eta(zone, k, j, i)
1530 + facteur(conte+2)*sol_hom3_eta(zone, k, j, i) ;
1532 mw.
set(zone, k, j, i) = sol_part_w(zone, k, j, i)
1533 + facteur(conte)*sol_hom1_w(zone, k, j, i)
1534 + facteur(conte+1)*sol_hom2_w(zone, k, j, i)
1535 + facteur(conte+2)*sol_hom3_w(zone, k, j, i) ;
1540 nr = mgrid.
get_nr(nzm1) ;
1541 for (
int i=0 ; i<nr ; i++) {
1542 mhrr.
set(nzm1, k, j, i) = sol_part_hrr(nzm1, k, j, i)
1543 + facteur(conte)*sol_hom1_hrr(nzm1, k, j, i)
1544 + facteur(conte+1)*sol_hom2_hrr(nzm1, k, j, i) ;
1546 meta.
set(nzm1, k, j, i) = sol_part_eta(nzm1, k, j, i)
1547 + facteur(conte)*sol_hom1_eta(nzm1, k, j, i)
1548 + facteur(conte+1)*sol_hom2_eta(nzm1, k, j, i) ;
1550 mw.
set(nzm1, k, j, i) = sol_part_w(nzm1, k, j, i)
1551 + facteur(conte)*sol_hom1_w(nzm1, k, j, i)
1552 + facteur(conte+1)*sol_hom2_w(nzm1, k, j, i) ;
1592 assert(mp_aff != 0x0) ;
1605 int nt = mgrid.
get_nt(0) ;
1606 int np = mgrid.
get_np(0) ;
1608 Scalar source = hh ;
1610 Scalar source_coq = hh ;
1612 source_coq.set_spectral_va().ylm() ;
1613 Base_val base = source.get_spectral_base() ;
1615 int lmax = base.give_lmax(mgrid, 0) + 1;
1622 Mtbl_cf sol_part_hrr(mgrid, base) ; sol_part_hrr.annule_hard() ;
1623 Mtbl_cf sol_part_eta(mgrid, base) ; sol_part_eta.annule_hard() ;
1624 Mtbl_cf sol_hom1_hrr(mgrid, base) ; sol_hom1_hrr.annule_hard() ;
1625 Mtbl_cf sol_hom1_eta(mgrid, base) ; sol_hom1_eta.annule_hard() ;
1626 Mtbl_cf sol_hom2_hrr(mgrid, base) ; sol_hom2_hrr.annule_hard() ;
1627 Mtbl_cf sol_hom2_eta(mgrid, base) ; sol_hom2_eta.annule_hard() ;
1629 bool need_calculation = true ;
1634 int l_q, m_q, base_r ;
1635 Itbl mat_done(lmax) ;
1641 int nr = mgrid.
get_nr(lz) ;
1642 double alpha = mp_aff->
get_alpha()[lz] ;
1643 Matrice ope(2*nr, 2*nr) ;
1644 if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1646 for (
int k=0 ; k<np+1 ; k++) {
1647 for (
int j=0 ; j<nt ; j++) {
1649 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
1650 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1651 if (need_calculation) {
1652 ope.set_etat_qcq() ;
1653 Diff_dsdx od(base_r, nr) ;
const Matrice& md = od.get_matrice() ;
1654 Diff_sx os(base_r, nr) ;
const Matrice& ms = os.get_matrice() ;
1656 for (
int lin=0; lin<nr; lin++)
1657 for (
int col=0; col<nr; col++)
1658 ope.set(lin,col) = md(lin,col) + 3*ms(lin,col) ;
1659 for (
int lin=0; lin<nr; lin++)
1660 for (
int col=0; col<nr; col++)
1661 ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
1662 for (
int lin=0; lin<nr; lin++)
1663 for (
int col=0; col<nr; col++)
1664 ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
1665 for (
int lin=0; lin<nr; lin++)
1666 for (
int col=0; col<nr; col++)
1667 ope.set(lin+nr,col+nr) = md(lin,col) + 3*ms(lin, col);
1670 for (
int col=0; col<2*nr; col++) {
1671 ope.set(nr-1, col) = 0 ;
1672 ope.set(2*nr-1, col) = 0 ;
1676 for (
int col=0; col<nr; col++) {
1677 ope.set(nr-1, col) = pari ;
1678 ope.set(2*nr-1, col+nr) = pari ;
1683 ope.set(nr-1, nr-1) = 1 ;
1684 ope.set(2*nr-1, 2*nr-1) = 1 ;
1688 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1689 Matrice* pope =
new Matrice(ope) ;
1691 mat_done.set(l_q) = 1 ;
1695 const Matrice& oper = (par_mat == 0x0 ? ope :
1698 sec.set_etat_qcq() ;
1699 for (
int lin=0; lin<nr; lin++)
1700 sec.set(lin) = (*source.get_spectral_va().c_cf)(lz, k, j, lin) ;
1701 for (
int lin=0; lin<nr; lin++)
1702 sec.set(nr+lin) = -0.5*(*source.get_spectral_va().c_cf)
1708 for (
int col=0; col<nr; col++) {
1710 (*source_coq.get_spectral_va().c_cf)(lz, k, j, col) ;
1713 sec.set(nr-1) = h0 / 3. ;
1715 sec.set(2*nr-1) = 0 ;
1716 Tbl sol = oper.inverse(sec) ;
1717 for (
int i=0; i<nr; i++) {
1718 sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1719 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1731 for (
int lz=1; lz<nz-1; lz++) {
1732 if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1733 int nr = mgrid.
get_nr(lz) ;
1736 assert(mgrid.
get_nt(lz) == nt) ;
1737 assert(mgrid.
get_np(lz) == np) ;
1738 double alpha = mp_aff->
get_alpha()[lz] ;
1739 double ech = mp_aff->
get_beta()[lz] / alpha ;
1740 Matrice ope(2*nr, 2*nr) ;
1742 for (
int k=0 ; k<np+1 ; k++) {
1743 for (
int j=0 ; j<nt ; j++) {
1745 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
1746 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1747 if (need_calculation) {
1748 ope.set_etat_qcq() ;
1749 Diff_xdsdx oxd(base_r, nr) ;
const Matrice& mxd = oxd.get_matrice() ;
1750 Diff_dsdx od(base_r, nr) ;
const Matrice& md = od.get_matrice() ;
1751 Diff_id oid(base_r, nr) ;
const Matrice& mid = oid.get_matrice() ;
1753 for (
int lin=0; lin<nr; lin++)
1754 for (
int col=0; col<nr; col++)
1755 ope.set(lin,col) = mxd(lin,col) + ech*md(lin,col)
1757 for (
int lin=0; lin<nr; lin++)
1758 for (
int col=0; col<nr; col++)
1759 ope.set(lin,col+nr) = -l_q*(l_q+1)*mid(lin,col) ;
1760 for (
int lin=0; lin<nr; lin++)
1761 for (
int col=0; col<nr; col++)
1762 ope.set(lin+nr,col) = -0.5*mid(lin,col) ;
1763 for (
int lin=0; lin<nr; lin++)
1764 for (
int col=0; col<nr; col++)
1765 ope.set(lin+nr,col+nr) = mxd(lin,col) + ech*md(lin,col)
1768 for (
int col=0; col<2*nr; col++) {
1769 ope.set(ind0+nr-1, col) = 0 ;
1770 ope.set(ind1+nr-1, col) = 0 ;
1772 ope.set(ind0+nr-1, ind0) = 1 ;
1773 ope.set(ind1+nr-1, ind1) = 1 ;
1776 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1777 Matrice* pope =
new Matrice(ope) ;
1779 mat_done.set(l_q) = 1 ;
1782 const Matrice& oper = (par_mat == 0x0 ? ope :
1785 sec.set_etat_qcq() ;
1786 for (
int lin=0; lin<nr; lin++)
1787 sec.set(lin) = (*source_coq.get_spectral_va().c_cf)
1789 for (
int lin=0; lin<nr; lin++)
1790 sec.set(nr+lin) = -0.5*(*source_coq.get_spectral_va().c_cf)
1792 sec.set(ind0+nr-1) = 0 ;
1793 sec.set(ind1+nr-1) = 0 ;
1794 Tbl sol = oper.inverse(sec) ;
1796 for (
int i=0; i<nr; i++) {
1797 sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1798 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1801 sec.set(ind0+nr-1) = 1 ;
1802 sol = oper.inverse(sec) ;
1803 for (
int i=0; i<nr; i++) {
1804 sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
1805 sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
1807 sec.set(ind0+nr-1) = 0 ;
1808 sec.set(ind1+nr-1) = 1 ;
1809 sol = oper.inverse(sec) ;
1810 for (
int i=0; i<nr; i++) {
1811 sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
1812 sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
1823 if (need_calculation && (par_mat != 0x0)) mat_done.annule_hard() ;
1824 int nr = mgrid.
get_nr(lz) ;
1827 assert(mgrid.
get_nt(lz) == nt) ;
1828 assert(mgrid.
get_np(lz) == np) ;
1829 double alpha = mp_aff->
get_alpha()[lz] ;
1830 Matrice ope(2*nr, 2*nr) ;
1832 for (
int k=0 ; k<np+1 ; k++) {
1833 for (
int j=0 ; j<nt ; j++) {
1835 base.give_quant_numbers(lz, k, j, m_q, l_q, base_r) ;
1836 if ( (nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1837 if (need_calculation) {
1838 ope.set_etat_qcq() ;
1839 Diff_dsdx od(base_r, nr) ;
const Matrice& md = od.get_matrice() ;
1840 Diff_sx os(base_r, nr) ;
const Matrice& ms = os.get_matrice() ;
1842 for (
int lin=0; lin<nr; lin++)
1843 for (
int col=0; col<nr; col++)
1844 ope.set(lin,col) = - md(lin,col) + 3*ms(lin,col) ;
1845 for (
int lin=0; lin<nr; lin++)
1846 for (
int col=0; col<nr; col++)
1847 ope.set(lin,col+nr) = -l_q*(l_q+1)*ms(lin,col) ;
1848 for (
int lin=0; lin<nr; lin++)
1849 for (
int col=0; col<nr; col++)
1850 ope.set(lin+nr,col) = -0.5*ms(lin,col) ;
1851 for (
int lin=0; lin<nr; lin++)
1852 for (
int col=0; col<nr; col++)
1853 ope.set(lin+nr,col+nr) = -md(lin,col) + 3*ms(lin, col) ;
1856 for (
int col=0; col<2*nr; col++) {
1857 ope.set(ind0+nr-2, col) = 0 ;
1858 ope.set(ind0+nr-1, col) = 0 ;
1859 ope.set(ind1+nr-2, col) = 0 ;
1860 ope.set(ind1+nr-1, col) = 0 ;
1862 for (
int col=0; col<nr; col++) {
1863 ope.set(ind0+nr-1, col+ind0) = 1 ;
1864 ope.set(ind1+nr-1, col+ind1) = 1 ;
1866 ope.set(ind0+nr-2, ind0+1) = 1 ;
1867 ope.set(ind1+nr-2, ind1+1) = 1 ;
1870 if ((par_mat != 0x0) && (mat_done(l_q) == 0)) {
1871 Matrice* pope =
new Matrice(ope) ;
1873 mat_done.set(l_q) = 1 ;
1876 const Matrice& oper = (par_mat == 0x0 ? ope :
1879 sec.set_etat_qcq() ;
1880 for (
int lin=0; lin<nr; lin++)
1881 sec.set(lin) = (*source.get_spectral_va().c_cf)
1883 for (
int lin=0; lin<nr; lin++)
1884 sec.set(nr+lin) = -0.5*(*source.get_spectral_va().c_cf)
1886 sec.set(ind0+nr-2) = 0 ;
1887 sec.set(ind0+nr-1) = 0 ;
1888 sec.set(ind1+nr-2) = 0 ;
1889 sec.set(ind1+nr-1) = 0 ;
1890 Tbl sol = oper.inverse(sec) ;
1891 for (
int i=0; i<nr; i++) {
1892 sol_part_hrr.set(lz, k, j, i) = sol(i) ;
1893 sol_part_eta.set(lz, k, j, i) = sol(i+nr) ;
1896 sec.set(ind0+nr-2) = 1 ;
1897 sol = oper.inverse(sec) ;
1898 for (
int i=0; i<nr; i++) {
1899 sol_hom1_hrr.set(lz, k, j, i) = sol(i) ;
1900 sol_hom1_eta.set(lz, k, j, i) = sol(i+nr) ;
1902 sec.set(ind0+nr-2) = 0 ;
1903 sec.set(ind1+nr-2) = 1 ;
1904 sol = oper.inverse(sec) ;
1905 for (
int i=0; i<nr; i++) {
1906 sol_hom2_hrr.set(lz, k, j, i) = sol(i) ;
1907 sol_hom2_eta.set(lz, k, j, i) = sol(i+nr) ;
1916 Mtbl_cf dhom1_hrr = sol_hom1_hrr ;
1917 Mtbl_cf dhom2_hrr = sol_hom2_hrr ;
1918 Mtbl_cf dpart_hrr = sol_part_hrr ;
1919 Mtbl_cf dhom1_eta = sol_hom1_eta ;
1920 Mtbl_cf dhom2_eta = sol_hom2_eta ;
1921 Mtbl_cf dpart_eta = sol_part_eta ;
1923 dhom1_hrr.dsdx() ; dhom1_eta.dsdx() ;
1924 dhom2_hrr.dsdx() ; dhom2_eta.dsdx() ;
1925 dpart_hrr.dsdx() ; dpart_eta.dsdx() ;
1929 int taille = 2*(nz-1) ;
1933 Tbl sec_membre(taille) ;
1934 Matrice systeme(taille, taille) ;
1935 int ligne ;
int colonne ;
1939 for (
int k=0 ; k<np+1 ; k++)
1940 for (
int j=0 ; j<nt ; j++) {
1941 base.give_quant_numbers(0, k, j, m_q, l_q, base_r) ;
1942 if ((nullite_plm(j, nt, k, np, base) == 1) && (l_q < 2)) {
1945 systeme.annule_hard() ;
1946 sec_membre.annule_hard() ;
1953 sec_membre.set(ligne) = 0.;
1956 sec_membre.set(ligne) = 0.;
1961 double alpha2 = mp_aff->
get_alpha()[1] ;
1968 systeme.set(ligne, colonne) = - (4. - double(l_q*(l_q+1.)))*sol_hom1_hrr.val_in_bound_jk(1,j,k) - 2.*dhom1_hrr.val_in_bound_jk(1,j,k)/alpha2;
1970 systeme.set(ligne, colonne +1) = - (4. - double(l_q*(l_q+1.)))*sol_hom2_hrr.val_in_bound_jk(1,j,k) - 2.*dhom2_hrr.val_in_bound_jk(1,j,k)/alpha2;
1972 sec_membre.set(ligne) = (4. - double(l_q*(l_q+1.)))*sol_part_hrr.val_in_bound_jk(1,j,k) + 2.*dpart_hrr.val_in_bound_jk(1,j,k)/alpha2 - blob;
1978 systeme.set(ligne, colonne) = - ( 2.*dhom1_eta.val_in_bound_jk(1,j,k)/alpha2 + (2.-double(l_q*(l_q+1.)))*sol_hom1_eta.val_in_bound_jk(1,j,k) + 2.*sol_hom1_hrr.val_in_bound_jk(1,j,k));
1980 systeme.set(ligne, colonne +1) = - ( 2.*dhom2_eta.val_in_bound_jk(1,j,k)/alpha2 + (2.-double(l_q*(l_q+1.)))*sol_hom2_eta.val_in_bound_jk(1,j,k) + 2.*sol_hom2_hrr.val_in_bound_jk(1,j,k));
1982 sec_membre.set(ligne) = 2.*dpart_eta.val_in_bound_jk(1,j,k)/alpha2 + (2.-double(l_q*(l_q+1.)))*sol_part_eta.val_in_bound_jk(1,j,k) + 2.*sol_part_hrr.val_in_bound_jk(1,j,k) - blob2;
1987 systeme.set(ligne, colonne) =
1988 sol_hom1_hrr.val_out_bound_jk(1, j, k) ;
1989 systeme.set(ligne, colonne+1) =
1990 sol_hom2_hrr.val_out_bound_jk(1, j, k) ;
1992 sec_membre.set(ligne) -= sol_part_hrr.val_out_bound_jk(1, j, k) ;
1995 systeme.set(ligne, colonne) =
1996 sol_hom1_eta.val_out_bound_jk(1, j, k) ;
1997 systeme.set(ligne, colonne+1) =
1998 sol_hom2_eta.val_out_bound_jk(1, j, k) ;
2000 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(1, j, k) ;
2006 for (
int zone=2 ; zone<nz-1 ; zone++) {
2007 nr = mgrid.
get_nr(zone) ;
2011 systeme.set(ligne, colonne) =
2012 - sol_hom1_hrr.val_in_bound_jk(zone, j, k) ;
2013 systeme.set(ligne, colonne+1) =
2014 - sol_hom2_hrr.val_in_bound_jk(zone, j, k) ;
2016 sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(zone, j, k) ;
2019 systeme.set(ligne, colonne) =
2020 - sol_hom1_eta.val_in_bound_jk(zone, j, k) ;
2021 systeme.set(ligne, colonne+1) =
2022 - sol_hom2_eta.val_in_bound_jk(zone, j, k) ;
2024 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(zone, j, k) ;
2028 systeme.set(ligne, colonne) =
2029 sol_hom1_hrr.val_out_bound_jk(zone, j, k) ;
2030 systeme.set(ligne, colonne+1) =
2031 sol_hom2_hrr.val_out_bound_jk(zone, j, k) ;
2033 sec_membre.set(ligne) -= sol_part_hrr.val_out_bound_jk(zone, j, k) ;
2036 systeme.set(ligne, colonne) =
2037 sol_hom1_eta.val_out_bound_jk(zone, j, k) ;
2038 systeme.set(ligne, colonne+1) =
2039 sol_hom2_eta.val_out_bound_jk(zone, j, k) ;
2041 sec_membre.set(ligne) -= sol_part_eta.val_out_bound_jk(zone, j, k) ;
2047 nr = mgrid.
get_nr(nz-1) ;
2051 systeme.set(ligne, colonne) =
2052 - sol_hom1_hrr.val_in_bound_jk(nz-1, j, k) ;
2053 systeme.set(ligne, colonne+1) =
2054 - sol_hom2_hrr.val_in_bound_jk(nz-1, j, k) ;
2056 sec_membre.set(ligne) += sol_part_hrr.val_in_bound_jk(nz-1, j, k) ;
2059 systeme.set(ligne, colonne) =
2060 - sol_hom1_eta.val_in_bound_jk(nz-1, j, k) ;
2061 systeme.set(ligne, colonne+1) =
2062 - sol_hom2_eta.val_in_bound_jk(nz-1, j, k) ;
2064 sec_membre.set(ligne) += sol_part_eta.val_in_bound_jk(nz-1, j, k) ;
2070 Tbl facteur = systeme.inverse(sec_membre) ;
2076 for (
int i=0 ; i<nr ; i++) {
2077 mhrr.set(0, k, j, i) = sol_part_hrr(0, k, j, i) ;
2078 meta.set(0, k, j, i) = sol_part_eta(0, k, j, i) ;
2080 for (
int zone=1 ; zone<nz-1 ; zone++) {
2081 nr = mgrid.
get_nr(zone) ;
2082 for (
int i=0 ; i<nr ; i++) {
2083 mhrr.set(zone, k, j, i) = sol_part_hrr(zone, k, j, i)
2084 + facteur(conte)*sol_hom1_hrr(zone, k, j, i)
2085 + facteur(conte+1)*sol_hom2_hrr(zone, k, j, i) ;
2087 meta.set(zone, k, j, i) = sol_part_eta(zone, k, j, i)
2088 + facteur(conte)*sol_hom1_eta(zone, k, j, i)
2089 + facteur(conte+1)*sol_hom2_eta(zone, k, j, i) ;
2093 nr = mgrid.
get_nr(nz-1) ;
2094 for (
int i=0 ; i<nr ; i++) {
2095 mhrr.set(nz-1, k, j, i) = sol_part_hrr(nz-1, k, j, i)
2096 + facteur(conte)*sol_hom1_hrr(nz-1, k, j, i)
2097 + facteur(conte+1)*sol_hom2_hrr(nz-1, k, j, i) ;
2099 meta.set(nz-1, k, j, i) = sol_part_eta(nz-1, k, j, i)
2100 + facteur(conte)*sol_hom1_eta(nz-1, k, j, i)
2101 + facteur(conte+1)*sol_hom2_eta(nz-1, k, j, i) ;
Bases of the spectral expansions.
void mult_x()
The basis is transformed as with a multiplication by .
int give_lmax(const Mg3d &mgrid, int lz) const
Returns the highest multipole for a given grid.
void give_quant_numbers(int, int, int, int &, int &, int &) const
Computes the various quantum numbers and 1d radial base.
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator Identity (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator division by (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Class for the elementary differential operator (see the base class Diff ).
virtual const Matrice & get_matrice() const
Returns the matrix associated with the operator.
Basic integer array class.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
int & set(int i)
Read/write of a particular element (index i ) (1D case)
void annule_hard()
Sets the Itbl to zero in a hard way.
const double * get_beta() const
Returns the pointer on the array beta.
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.
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
double & set(int j, int i)
Read/write of a particuliar element.
Tbl inverse(const Tbl &sec_membre) const
Solves the linear system represented by the matrix.
void annule_hard()
Sets the logical state to ETATQCQ (undefined state).
void set_lu() const
Calculate the LU-representation, assuming the band-storage has been done.
int get_type_t() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the equatorial pl...
int get_np(int l) const
Returns the number of points in the azimuthal direction ( ) in domain no. l.
int get_type_p() const
Returns the type of sampling in the direction: SYM : : symmetry with respect to the transformatio...
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.
Coefficients storage for the multi-domain spectral method.
void annule_hard()
Sets the Mtbl_cf to zero in a hard way.
double val_out_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
double val_in_bound_jk(int l, int j, int k) const
Computes the angular coefficient of index j,k of the field represented by *this at by means of the s...
Tbl & set(int l)
Read/write of the Tbl containing the coefficients in a given domain.
int get_n_itbl_mod() const
Returns the number of modifiable Itbl 's addresses in the list.
int get_n_matrice_mod() const
Returns the number of modifiable Matrice 's addresses in the list.
const int & get_int(int position=0) const
Returns the reference of a int stored in the list.
int get_n_tbl_mod() const
Returns the number of modifiable Tbl 's addresses in the list.
Itbl & get_itbl_mod(int position=0) const
Returns the reference of a stored modifiable Itbl .
int get_n_int_mod() const
Returns the number of modifiable int 's addresses in the list.
Tbl & get_tbl_mod(int position=0) const
Returns the reference of a modifiable Tbl stored in the list.
void clean_all()
Deletes all the objects stored as modifiables, i.e.
void add_matrice_mod(Matrice &ti, int position=0)
Adds the address of a new modifiable Matrice to the list.
Matrice & get_matrice_mod(int position=0) const
Returns the reference of a modifiable Matrice stored in the list.
void add_int_mod(int &n, int position=0)
Adds the address of a new modifiable int to the list.
void add_tbl_mod(Tbl &ti, int position=0)
Adds the address of a new modifiable Tbl to the list.
void add_itbl_mod(Itbl &ti, int position=0)
Adds the address of a new modifiable Itbl to the list.
int get_n_int() const
Returns the number of stored int 's addresses.
int & get_int_mod(int position=0) const
Returns the reference of a modifiable int stored in the list.
Tensor field of valence 0 (or component of a tensorial field).
void annule_l(int l_min, int l_max, bool ylm_output=false)
Sets all the multipolar components between l_min and l_max to zero.
const Base_val & get_spectral_base() const
Returns the spectral bases of the Valeur va
void div_r_dzpuis(int ced_mult_r)
Division by r everywhere but with the output flag dzpuis set to ced_mult_r .
virtual void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
bool check_dzpuis(int dzi) const
Returns false if the last domain is compactified and *this is not zero in this domain and dzpuis is n...
const Scalar & dsdr() const
Returns of *this .
void annule_hard()
Sets the Scalar to zero in a hard way.
int get_etat() const
Returns the logical state ETATNONDEF (undefined), ETATZERO (null) or ETATQCQ (ordinary).
const Valeur & get_spectral_va() const
Returns va (read only version)
Valeur & set_spectral_va()
Returns va (read/write version)
void mult_r_dzpuis(int ced_mult_r)
Multiplication by r everywhere but with the output flag dzpuis set to ced_mult_r .
void mult_r()
Multiplication by r everywhere; dzpuis is not changed.
void set_spectral_base(const Base_val &)
Sets the spectral bases of the Valeur va
void sol_Dirac_A2(const Scalar &aaa, Scalar &tilde_mu, Scalar &x_new, Scalar bound_mu, const Param *par_bc)
Same resolution as sol_Dirac_Abound, but here the boundary conditions are the degenerate elliptic con...
void sol_Dirac_BC3(const Scalar &bb, const Scalar &hh, Scalar &hrr, Scalar &tilde_eta, Scalar &ww, Scalar bound_hrr, Scalar bound_eta, Param *par_bc, Param *par_mat)
Same resolution as sol_Dirac_Abound, but here the boundary conditions are the degenerate elliptic con...
void annule_hard()
Sets the Tbl to zero in a hard way.
void set_etat_zero()
Sets the logical state to ETATZERO (zero).
double & set(int i)
Read/write of a particular element (index i) (1D case)
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
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 * c
Values of the function at the points of the multi-grid
Mtbl_cf * c_cf
Coefficients of the spectral expansion of the function.
void ylm_i()
Inverse of ylm()
#define R_CHEBP
base de Cheb. paire (rare) seulement
const Map *const mp
Mapping on which the numerical values at the grid points are defined.
void annule_domain(int l)
Sets the Tensor to zero in a given domain.