LORENE
prolonge_c1.C
1/*
2 * Small utility for determining the surface for 2-fluid stars.
3 *
4 */
5
6/*
7 * Copyright (c) 2002 Jerome Novak
8 *
9 * This file is part of LORENE.
10 *
11 * LORENE is free software; you can redistribute it and/or modify
12 * it under the terms of the GNU General Public License version 2
13 * as published by the Free Software Foundation.
14 *
15 * LORENE is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the GNU General Public License
21 * along with LORENE; if not, write to the Free Software
22 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 *
24 */
25
26char prolonge_c1_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Utilities/prolonge_c1.C,v 1.4 2016/09/19 15:26:23 j_novak Exp $" ;
27
28#include "cmp.h"
29
30namespace Lorene {
31Cmp prolonge_c1(const Cmp& uu, const int nzet) {
32
33 const Map_radial* mpi = dynamic_cast<const Map_radial*>( uu.get_mp() ) ;
34
35 if (mpi == 0x0) {
36 cout <<
37 "prolonge_c1 : The mapping does not belong to the class Map_radial !"
38 << endl ;
39 abort() ;
40 }
41
42 const Map_radial& mp = *mpi ;
43
44 const Mg3d& mg = *(mp.get_mg()) ;
45
46 int nz = mg.get_nzone() ;
47 assert((nzet > 0)&&(nzet<nz)) ;
48
49 Cmp resu(mp) ;
50 resu.allocate_all() ;
51 double nbc = uu(0, 0, 0, 0) ;
52 double resu_ext = - 0.2 * nbc ;
53 const Coord& rr = mp.r ;
54 double rout = (+rr)(nzet,0,0,mg.get_nr(nzet)-1) ;
55 double rin = 0 ; double resu1 = 0 ; double dresu1 = 0 ;
56 double a = 0 ; double b = 0 ;
57 for (int k=0; k<mg.get_np(0); k++)
58 for (int j=0; j<mg.get_nt(0); j++) {
59 double ir = 0 ;
60 double nm2 = nbc ;
61 double nm1 = nbc ;
62 double nb0 = nbc ;
63 double np1 = nbc ;
64 double rm2 = 0 ;
65 double rm1 = 0 ;
66 double r0 = 0 ;
67 double rp1 = 0 ;
68 bool dedans = true ;
69 for (int lz=0; lz <= nzet; lz++) {
70 for (int i=1; i<mg.get_nr(lz); i++) {
71 ir++ ;
72 rm2 = rm1 ;
73 rm1 = r0 ;
74 r0 = rp1 ;
75 rp1 = (+rr)(lz,k,j,i) ;
76 nm2 = nm1 ;
77 nm1 = nb0 ;
78 nb0 = np1 ;
79 np1 = uu(lz,k,j,i) ;
80 if ((np1<=0.) && dedans) {
81 if (ir<2) {
82 cout << "Problem prolonge_c1!" << endl ;
83 abort() ;
84 }
85 resu1 = nm1 * (r0-rp1)*(rm2-rp1)
86 / ((r0-rm1)*(rm2-rm1))
87 + nm2 * (r0-rp1)*(rm1-rp1) / ((r0-rm2)*(rm1-rm2))
88 + nb0 * (rm1-rp1)*(rm2-rp1) / ((rm1-r0)*(rm2-r0)) ;
89 resu.set(lz,k,j,i) = resu1 ;
90 dresu1 = nm1 * ((rp1-r0) + (rp1-rm2))
91 / ((r0-rm1)*(rm2-rm1))
92 + nm2 * ((rp1-r0) + (rp1-rm1)) / ((r0-rm2)*(rm1-rm2))
93 + nb0 * ((rp1-rm1) + (rp1-rm2)) / ((rm1-r0)*(rm2-r0)) ;
94 a = (dresu1 - 2*(resu_ext - resu1)/(rout - rp1))/
95 ((rout-rp1)*(rout-rp1)) ;
96 b = 0.5*(-dresu1/(rout-rp1) - a*(rout+rp1)) ;
97 rin = rp1 ;
98 dedans = false ;
99 }
100 else {
101 resu.set(lz,k,j,i) = (dedans ? np1 :
102 resu1*(rout-rp1)/(rout-rin) + resu_ext*(rin-rp1)/(rin-rout)
103 +(rp1-rin)*(rp1-rout)*(a*rp1+b) );
104 }
105 }
106 resu.set(lz,k,j,0) = (lz==0 ? nbc :
107 resu(lz-1, k, j, mg.get_nr(lz-1)-1 ) ) ;
108 }
109 }
110 Cmp resu2(mp) ;
111 resu2 = resu_ext ;
112 resu2.annule(0,nzet) ;
113 resu.annule(nzet+1, nz-1) ;
114 resu += resu2 ;
115 resu.std_base_scal() ;
116 return resu ;
117}
118}
Lorene prototypes.
Definition app_hor.h:64