Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Open sidebar
Jean-Eric Campagne
LagSHT
Commits
f0d1bbd2
Commit
f0d1bbd2
authored
Apr 09, 2015
by
Jean-Eric Campagne
Browse files
(JEC) 9/4/15 delete Sophya old stuff
parent
f0acd5a6
Changes
9
Expand all
Hide whitespace changes
Inline
Side-by-side
Showing
9 changed files
with
0 additions
and
2466 deletions
+0
-2466
Sophya/Makefile
Sophya/Makefile
+0
-73
Sophya/extrafunctions.h
Sophya/extrafunctions.h
+0
-33
Sophya/lagSphericTransform.cc
Sophya/lagSphericTransform.cc
+0
-321
Sophya/lagSphericTransform.h
Sophya/lagSphericTransform.h
+0
-142
Sophya/laguerreBuilder.cc
Sophya/laguerreBuilder.cc
+0
-215
Sophya/laguerreBuilder.h
Sophya/laguerreBuilder.h
+0
-119
Sophya/laguerreTransform.cc
Sophya/laguerreTransform.cc
+0
-448
Sophya/laguerreTransform.h
Sophya/laguerreTransform.h
+0
-125
Sophya/tstlag.cc
Sophya/tstlag.cc
+0
-990
No files found.
Sophya/Makefile
deleted
100644 → 0
View file @
f0acd5a6
include
$(SOPHYABASE)/include/sophyamake.inc
# ===== Libsharp =====
SHARPDIR
=
/Users/campagne/Travail/kits/libsharp-code/auto/
SHARPLIB
=
$(SHARPDIR)
/lib
SHARPINC
=
$(SHARPDIR)
/include
SHARPLIBN
=
-lsharp
-lc_utils
-lfftpack
# ===== NTL
NTLLIB
=
/opt/local/lib/libntl.a
CPPFLAGS
+=
-I
$(SHARPINC)
-I
/opt/local/include
LDFLAGS
+=
-L
$(SHARPLIB)
$(SHARPLIBN)
$(NTLLIB)
-lm
OBJ
=
./Objs/
EXE
=
./Objs/
# Define our target list
.PHONY
:
default
default
:
tstlag
.PHONY
:
all
all
:
tstlag
.PHONY
:
tidy
tidy
:
rm
*
~
.PHONY
:
clean
clean
:
rm
-f
./Objs/
*
#C++ common Objects
CXXOBJ
=
$(OBJ)
laguerreBuilder.o
\
$(OBJ)
laguerreTransform.o
\
$(OBJ)
lagSphericTransform.o
#C++ common Headers
CXXHDR
=
laguerreBuilder.h
\
laguerreTransform.h
\
lagSphericTransform.h
\
extrafunctions.h
CXXFLAGS
=
-Wa
,-q
-march
=
native
-fopenmp
-fno-common
-O3
CXXCOMPILE
=
$(CXX)
$(CPPFLAGS)
$(CXXFLAGS)
-c
#LDFLAGS += $(OPENMP) -lm
#C++ rule for compiling
$(OBJ)%.o
:
%.cc $(CXXHDR)
echo
"compile...
$@
depend
$^
"
$(CXXCOMPILE)
-c
$<
-o
$@
######################
.PHONY
:
tstlag
tstlag
:
$(EXE)tstlag
echo
'---- tstlag made'
$(EXE)tstlag
:
$(OBJ)tstlag.o $(CXXOBJ)
echo
"Link..."
$(CXXLINK)
-o
$@
$(OBJ)
tstlag.o
$(CXXOBJ)
$(LDFLAGS)
$(SOPHYAEXTSLBLIST)
$(OBJ)tstlag.o
:
tstlag.cc $(CXXHDR)
echo
"compile...
$@
"
$(CXXCOMPILE)
$(FFLAGS)
-c
$<
-o
$@
Sophya/extrafunctions.h
deleted
100644 → 0
View file @
f0acd5a6
#ifndef EXTRAFUNCTIONS_SEEN
#define EXTRAFUNCTIONS_SEEN
#ifdef SO_NTLXDBLE
#include "NTL/xdouble.h"
/* double : mantisse, long: exponent (NTL-library) */
using
namespace
NTL
;
#endif
namespace
SOPHYA
{
//Utility fonctions to unify the code writing between "long double" (r_16) and NTL/xdouble
#ifdef SO_NTLXDBLE
inline
xdouble
exp
(
const
xdouble
&
x
)
{
return
xexp
(
to_double
(
x
));
}
#endif
inline
long
double
exp
(
long
double
x
)
{
return
expl
(
x
);
}
inline
long
double
fabs
(
long
double
x
)
{
return
fabsl
(
x
);
}
inline
double
to_double
(
long
double
x
)
{
return
double
(
x
);
}
}
//End namespace
#endif
Sophya/lagSphericTransform.cc
deleted
100644 → 0
View file @
f0acd5a6
#include <stdlib.h>
#include <iostream>
//Sophya
#include "sopnamsp.h"
#include "machdefs.h"
#include "pexceptions.h"
#define LAGSPHERICTRANSFORM_CC_BFILE // avoid extern template declarations
#include "lagSphericTransform.h"
#define DEBUG 0
namespace
SOPHYA
{
void
BaseLagSphTransform
::
SetThetaPhiMap
(
string
choice
){
//Todo: on pourrait choisir ici son type de geoemetrie 2D: ici ECP = Fejer's 1 sum rule
//libsharp defini Lmax = L-1 avec une somme l:0...Lmax inclus
//dans sharp_testsuite.c Fejer's 1st rule
//nrings = 2Lmax+1
//npixperring = 2Mmax+1
//et pour les Alm
//set_trinagular_alm (Lmax,Mmax).
Nrings_
=
2
*
L_
-
1
;
Nphi_
=
2
*
M_
-
1
;
sharpjob_
.
set_ECP_geometry
(
Nrings_
,
Nphi_
);
sharpjob_
.
set_triangular_alm_info
(
L_
-
1
,
M_
-
1
);
Npix_
=
get_npix
(
sharpjob_
.
get_geom_info
());
Nalm_
=
get_nalms
(
sharpjob_
.
get_alm_info
());
geomDone_
=
true
;
cout
<<
"SetThetaPhiMap: Nring, Nphi, Npix, Nalm "
<<
Nrings_
<<
", "
<<
Nphi_
<<
", "
<<
Npix_
<<
", "
<<
Nalm_
<<
endl
;
}
//SetThetaPhiMap
template
<
class
T
>
void
LaguerreSphericalTransform
<
T
>::
Synthesis
(
const
NDataBlock
<
complex
<
T
>
>&
flmn
,
NDataBlock
<
T
>&
fijk
)
{
int
Ntot3DPix
=
Npix_
*
N_
;
int
NtotCoeff
=
Nalm_
*
N_
;
if
(
flmn
.
Size
()
!=
NtotCoeff
)
throw
SzMismatchError
(
"SphLagTransform::Synthesis size error"
);
fijk
.
ReSize
(
Ntot3DPix
);
//Perform the Inverse Laguerre Transform first
NDataBlock
<
complex
<
T
>
>
xflmk
(
NtotCoeff
);
//T version of the flmk
#if DEBUG >= 1
cout
<<
"SphLagTransform::Synthesis type <T> I/O .... Go..."
<<
endl
;
cout
<<
"Start Inverse Laguerre...."
<<
endl
;
#endif
lagTrans_
.
MultiSynthesis
(
flmn
,
xflmk
,
Nalm_
);
#if DEBUG >= 1
cout
<<
"End Inverse Laguerre...."
<<
endl
;
cout
<<
"Start Inverse H.Spherical Transform "
<<
endl
;
#endif
int
off7n
,
off7k
;
//perform the synthesis per shell
NDataBlock
<
complex
<
r_8
>
>
flmk
(
Nalm_
);
// r_8 version of the flmk for libsharp
NDataBlock
<
T
>
nodes
(
lagTrans_
.
GetNodes
());
for
(
int
n
=
0
;
n
<
N_
;
n
++
){
off7n
=
n
*
Nalm_
;
off7k
=
n
*
Npix_
;
T
scale
=
exp
(
-
nodes
(
n
)
/
2.
);
//<-------- use a scale to call sharp code
#if DEBUG >= 1
cout
<<
"Synthesis... at shell num: "
<<
n
<<
" scale = "
<<
scale
<<
endl
;
#endif
//#pragma omp parallel for schedule(dynamic,1)
for
(
int
i
=
0
;
i
<
Nalm_
;
i
++
){
xflmk
(
off7n
+
i
)
*=
scale
;
flmk
(
i
)
=
complex
<
r_8
>
(
to_double
(
xflmk
(
off7n
+
i
).
real
()),
to_double
(
xflmk
(
off7n
+
i
).
imag
()));
#if DEBUG >= 2
cout
<<
" Syntheis... flmk("
<<
i
<<
"): "
<<
xflmk
(
off7n
+
i
)
<<
", "
<<
flmk
(
i
)
<<
endl
;
#endif
}
NDataBlock
<
r_8
>
rfijk
(
Npix_
);
//version for libsharp : todo see if declaration outside the loops is better
sharpjob_
.
alm2map
(
&
(
flmk
.
Data
()[
0
].
real
()),
&
rfijk
.
Data
()[
0
],
false
);
//the false is to forbidd accumulation
//#pragma omp parallel for schedule(dynamic,1)
for
(
int
i
=
0
;
i
<
Npix_
;
i
++
){
fijk
(
off7k
+
i
)
=
rfijk
(
i
);
//store & switch r_8 -> T
fijk
(
off7k
+
i
)
/=
scale
;
//rescaling
#if DEBUG >= 2
cout
<<
" Syntheis... rescaled fijk("
<<
i
<<
"): "
<<
fijk
(
off7k
+
i
)
<<
" rfijk: "
<<
rfijk
(
i
)
<<
endl
;
#endif
}
}
//loop on shell
#if DEBUG >= 1
cout
<<
"SphLagTransform::Synthesis.... Done"
<<
endl
;
#endif
}
//Synthesis
template
<
class
T
>
void
LaguerreSphericalTransform
<
T
>::
Analysis
(
const
NDataBlock
<
T
>&
fijk
,
NDataBlock
<
complex
<
T
>
>&
flmn
)
{
int
Ntot3DPix
=
Npix_
*
N_
;
int
NtotCoeff
=
Nalm_
*
N_
;
if
(
fijk
.
Size
()
!=
Ntot3DPix
)
throw
SzMismatchError
(
"SphLagTransform::Analysis size error"
);
flmn
.
ReSize
(
NtotCoeff
);
#if DEBUG >= 1
cout
<<
"SphLagTransform::Analysis xdbl.... Go..."
<<
endl
;
#endif
//Perform the Spherical Hamonic analysis first
NDataBlock
<
complex
<
T
>
>
xflmk
(
NtotCoeff
);
NDataBlock
<
T
>
nodes
(
lagTrans_
.
GetNodes
());
NDataBlock
<
complex
<
r_8
>
>
flmk
(
Nalm_
);
//used by lisharp
NDataBlock
<
r_8
>
rfijk
(
Npix_
);
//used by libsharp
int
off7n
,
off7k
;
//Ylm analysis per shell
for
(
int
n
=
0
;
n
<
N_
;
n
++
){
off7n
=
n
*
Nalm_
;
off7k
=
n
*
Npix_
;
#if DEBUG >= 1
cout
<<
"Start Harmonic Spheric... at n = "
<<
n
<<
endl
;
#endif
//rescalling des pixels values
T
scale
=
exp
(
nodes
(
n
)
/
2.
);
//#pragma omp parallel for schedule(dynamic,1)
for
(
int
i
=
0
;
i
<
Npix_
;
i
++
){
rfijk
(
i
)
=
to_double
(
fijk
(
off7k
+
i
)
*
scale
);
//store & switch T -> r_8 with a scale
#if DEBUG >= 2
cout
<<
"rfijk("
<<
i
<<
"): "
<<
rfijk
(
i
)
<<
", "
<<
fijk
(
off7k
+
i
)
<<
endl
;
#endif
}
sharpjob_
.
map2alm
(
&
(
rfijk
.
Data
()[
0
]),
&
(
flmk
.
Data
()[
0
].
real
()),
false
);
//the false is to forbidd accumulation
//#pragma omp parallel for schedule(dynamic,1)
for
(
int
i
=
0
;
i
<
Nalm_
;
i
++
){
xflmk
(
i
+
off7n
)
=
flmk
(
i
);
xflmk
(
i
+
off7n
)
/=
scale
;
// rescaling
#if DEBUG >= 2
cout
<<
"flmk("
<<
i
<<
"): "
<<
xflmk
(
i
+
off7n
)
<<
", "
<<
flmk
(
i
+
off7n
)
<<
endl
;
#endif
}
}
//loop on shell
#if DEBUG >= 1
cout
<<
"End Harmonic Spheric..."
<<
endl
;
cout
<<
"Start Laguerre..."
<<
endl
;
#endif
//Laguerre Transform
lagTrans_
.
MultiAnalysis
(
xflmk
,
flmn
,
Nalm_
);
#if DEBUG >= 1
cout
<<
"End Laguerre..."
<<
endl
;
#endif
#if DEBUG >= 2
for
(
int
i
=
0
;
i
<
NtotCoeff
;
i
++
)
{
cout
<<
"flmn("
<<
i
<<
"): "
<<
flmn
(
i
)
<<
endl
;
}
#endif
#if DEBUG >= 1
cout
<<
"SphLagTransform::Analysis.... Done"
<<
endl
;
#endif
}
//Analysis
///////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////
#ifdef __CXX_PRAGMA_TEMPLATES__
#ifdef SO_LDBLE128
#pragma define_template LaguerreSphericalTransform<r_16>
#endif
#ifdef SO_NTLXDBLE
#pragma define_template LaguerreSphericalTransform<xdouble>
#endif
#endif
#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
#ifdef SO_LDBLE128
template
class
LaguerreSphericalTransform
<
r_16
>;
#endif
#ifdef SO_NTLXDBLE
template
class
LaguerreSphericalTransform
<
xdouble
>;
#endif
#endif
template
<
>
LaguerreSphericalTransform
<
r_8
>::
LaguerreSphericalTransform
(
int
Lmax
,
int
Mmax
,
int
Nmax
,
r_8
R
,
int
alpha
)
:
BaseLagSphTransform
(
Lmax
,
Mmax
,
Nmax
,
alpha
),
R_
(
R
),
lagTrans_
(
Nmax
,
R
,
alpha
)
{
}
//Ctor
template
<
>
LaguerreSphericalTransform
<
r_8
>::~
LaguerreSphericalTransform
()
{}
template
<
>
void
LaguerreSphericalTransform
<
r_8
>::
Synthesis
(
const
NDataBlock
<
complex
<
r_8
>
>&
flmn
,
NDataBlock
<
r_8
>&
fijk
)
{
int
Ntot3DPix
=
Npix_
*
N_
;
int
NtotCoeff
=
Nalm_
*
N_
;
if
(
flmn
.
Size
()
!=
NtotCoeff
)
throw
SzMismatchError
(
"SphLagTransform::Synthesis size error"
);
fijk
.
ReSize
(
Ntot3DPix
);
//Perform the Inverse Laguerre Transform first
NDataBlock
<
complex
<
r_8
>
>
flmk
(
NtotCoeff
);
#if DEBUG >= 1
cout
<<
"SphLagTransform::Synthesis type <r_8> I/O .... Go..."
<<
endl
;
cout
<<
"Start Inverse Laguerre...."
<<
endl
;
#endif
lagTrans_
.
MultiSynthesis
(
flmn
,
flmk
,
Nalm_
);
#if DEBUG >= 1
cout
<<
"End Inverse Laguerre...."
<<
endl
;
cout
<<
"Start Inverse H.Spherical Transform "
<<
endl
;
#endif
int
off7n
,
off7k
;
//perform the synthesis per shell
for
(
int
n
=
0
;
n
<
N_
;
n
++
){
off7n
=
n
*
Nalm_
;
off7k
=
n
*
Npix_
;
#if DEBUG >= 2
cout
<<
"Synthesis... at shell num: "
<<
n
<<
endl
;
#endif
sharpjob_
.
alm2map
(
&
(
flmk
.
Data
()[
off7n
].
real
()),
&
fijk
.
Data
()[
off7k
],
false
);
}
//loop on shell
#if DEBUG >= 1
cout
<<
"SphLagTransform::Synthesis.... Done"
<<
endl
;
#endif
}
//Synthesis
template
<
>
void
LaguerreSphericalTransform
<
r_8
>::
Analysis
(
const
NDataBlock
<
r_8
>&
fijk
,
NDataBlock
<
complex
<
r_8
>
>&
flmn
)
{
int
Ntot3DPix
=
Npix_
*
N_
;
int
NtotCoeff
=
Nalm_
*
N_
;
if
(
fijk
.
Size
()
!=
Ntot3DPix
)
throw
SzMismatchError
(
"SphLagTransform::Analysis size error"
);
flmn
.
ReSize
(
NtotCoeff
);
#if DEBUG >= 1
cout
<<
"SphLagTransform::Analysis xdbl.... Go..."
<<
endl
;
#endif
//Perform the Spherical Hamonic analysis first
NDataBlock
<
complex
<
r_8
>
>
flmk
(
NtotCoeff
);
int
off7n
,
off7k
;
//Ylm analysis per shell
for
(
int
n
=
0
;
n
<
N_
;
n
++
){
off7n
=
n
*
Nalm_
;
off7k
=
n
*
Npix_
;
#if DEBUG >= 1
cout
<<
"Start Harmonic Spheric... at n = "
<<
n
<<
endl
;
#endif
// sharpjob_.map2alm(&(fijk.Data()[0]),&(flmk.Data()[0].real()),false); //the false is to forbidd accumulation
sharpjob_
.
map2alm
(
&
(
fijk
.
Data
()[
off7k
]),
&
(
flmk
.
Data
()[
off7n
].
real
()),
false
);
//the false is to forbidd accumulation
}
//loop on shell
#if DEBUG >= 1
cout
<<
"End Harmonic Spheric..."
<<
endl
;
cout
<<
"Start Laguerre..."
<<
endl
;
#endif
//Laguerre Transform
lagTrans_
.
MultiAnalysis
(
flmk
,
flmn
,
Nalm_
);
#if DEBUG >= 1
cout
<<
"End Laguerre..."
<<
endl
;
#endif
#if DEBUG >= 2
for
(
int
i
=
0
;
i
<
NtotCoeff
;
i
++
)
{
cout
<<
"flmn("
<<
i
<<
"): "
<<
flmn
(
i
)
<<
endl
;
}
#endif
#if DEBUG >= 1
cout
<<
"SphLagTransform::Analysis.... Done"
<<
endl
;
#endif
}
//Analysis
}
// Fin de namespace
Sophya/lagSphericTransform.h
deleted
100644 → 0
View file @
f0acd5a6
#ifndef LAGSPHERICTRANSFORM_SEEN
#define LAGSPHERICTRANSFORM_SEEN
//Sophya
#include "machdefs.h"
#include "ndatablock.h"
//Sharp
#include "sharp_lowlevel.h"
#include "sharp_geomhelpers.h"
#include "sharp_almhelpers.h"
#include "sharp_cxx.h"
// extern "C" {
// #include <complex.h>
// }
#include "laguerreTransform.h"
#include "extrafunctions.h"
namespace
SOPHYA
{
////////////////////////////////////////////////////////////////
//// --------- Class Spherical Harmonic Laguerre Transform -- //
////////////////////////////////////////////////////////////////
class
BaseLagSphTransform
{
public:
//! Creator
BaseLagSphTransform
(
int
Lmax
,
int
Mmax
,
int
Nmax
,
int
alpha
=
2
)
:
L_
(
Lmax
),
M_
(
Mmax
),
N_
(
Nmax
),
alpha_
(
alpha
),
Nrings_
(
0
),
Nphi_
(
0
),
Npix_
(
0
),
Nalm_
(
0
),
geomDone_
(
false
)
{
if
(
Mmax
<
0
)
M_
=
Lmax
;
}
//! Destructor
virtual
~
BaseLagSphTransform
()
{}
//! Number of colatitude rings of the 2D-sphere representation
int
Nrings
()
{
return
Nrings_
;
}
//! Number of colatitude rings of the 2D-sphere representation
int
NTheta
()
{
return
Nrings
();
}
//! Number of phi values (pixels) per colatitude rings
int
NPhi
()
{
return
Nphi_
;
}
//! Total number of 2D pixels
int
Npix
()
{
return
Npix_
;
}
//! Total number of Alm coefficients
int
Nalm
()
{
return
Nalm_
;
}
//! Choose the 2D map representation & initialization of libsharp
void
SetThetaPhiMap
(
string
choice
=
""
);
protected:
//from sharp test suite
static
ptrdiff_t
get_npix
(
const
sharp_geom_info
*
ginfo
)
{
ptrdiff_t
res
=
0
;
for
(
int
i
=
0
;
i
<
ginfo
->
npairs
;
++
i
){
res
+=
ginfo
->
pair
[
i
].
r1
.
nph
;
if
(
ginfo
->
pair
[
i
].
r2
.
nph
>
0
)
res
+=
ginfo
->
pair
[
i
].
r2
.
nph
;
}
return
res
;
}
//from sharp test suite
static
ptrdiff_t
get_nalms
(
const
sharp_alm_info
*
ainfo
)
{
ptrdiff_t
res
=
0
;
for
(
int
i
=
0
;
i
<
ainfo
->
nm
;
++
i
)
res
+=
ainfo
->
lmax
-
ainfo
->
mval
[
i
]
+
1
;
return
res
;
}
protected:
int
L_
;
//!< Spherical harmonic limit such that l:0,...,L-1 m:0,...,l
int
M_
;
//!< maximum value of m
int
N_
;
//!< Radial band limit such that n:0,...,N-1
int
alpha_
;
//!< parameter of the generalized Laguerre Polynomials
int
Nrings_
;
//!< Number of colatitude rings of the 2D-sphere representation
int
Nphi_
;
//!< Number of phi values (pixels) per colatitude rings
int
Npix_
;
//!< Total number of 2D pixels
int
Nalm_
;
//!< Total number of Alm coefficients
bool
geomDone_
;
//!< Set to true once the geometry is properly set (See SetThetaPhiMap(string) )
sharp_cxxjob
<
r_8
>
sharpjob_
;
//!< libsharp job // voir si r_8 -> r_16 marcherait !!!
};
template
<
class
T
>
class
LaguerreSphericalTransform
:
public
BaseLagSphTransform
{
public:
//! Constructor
LaguerreSphericalTransform
(
int
Lmax
,
int
Mmax
,
int
Nmax
,
r_8
R
,
int
alpha
=
2
)
:
BaseLagSphTransform
(
Lmax
,
Mmax
,
Nmax
,
alpha
),
R_
(
R
),
lagTrans_
(
Nmax
,(
T
)
R
,
alpha
)
{
}
//Ctor
//! Destructor
virtual
~
LaguerreSphericalTransform
()
{}
//! Synthesis
/*! \brief Coeffs -> Pixels with Input/Output using T floating representation
\param flmn: 3D complex spherical-laguerre coefficients
\param fijk: 3D real function values
*/
void
Synthesis
(
const
NDataBlock
<
complex
<
T
>
>&
flmn
,
NDataBlock
<
T
>&
fijk
);
//! Analysis
/*! \brief Pixels -> Coeffs with Input/Output using T floating representation
\param fijk: 3D real function values
\param flmn: 3D complex spherical-laguerre coefficients
*/
void
Analysis
(
const
NDataBlock
<
T
>&
fijk
,
NDataBlock
<
complex
<
T
>
>&
flmn
);
protected:
r_8
R_
;
//!< Radial maximal value (not yet used)
LaguerreTransform
<
T
>
lagTrans_
;
//!< the Laguerre Transform
};
//class
//--------- extern template declarations (if needed) -----------
#if defined ( NEED_EXT_DECL_TEMP ) && !defined( LAGSPHERICTRANSFORM_CC_BFILE )
extern
template
class
LaguerreSphericalTransform
<
r_8
>;
#ifdef SO_LDBLE128
extern
template
class
LaguerreSphericalTransform
<
r_16
>;
#endif
#ifdef SO_NTLXDBLE
extern
template
class
LaguerreSphericalTransform
<
xdouble
>;
#endif
#endif // Fin de if defined ( NEED_EXT_DECL_TEMP )
}
//fin du namesape
#endif // LAGSPHERICTRANSFORM
Sophya/laguerreBuilder.cc
deleted
100644 → 0
View file @
f0acd5a6
#include <stdlib.h>
#include <iostream>
//Sophya
#include "sopnamsp.h"
#include "machdefs.h"
#include "pexceptions.h"
#define LAGUERREBUILDER_CC_BFILE // avoid extern template declarations
#include "laguerreBuilder.h"
#define DEBUG 1
namespace
SOPHYA
{
template
<
class
T
>
void
LaguerrePolyQuad
<
T
>::
QuadWeightNodes
(
NDataBlock
<
T
>&
nodes
,
NDataBlock
<
T
>&
weights
){
nodes
.
ReSize
(
N_
);
weights
.
ReSize
(
N_
);
T
x
;
T
tN
=
(
T
)
N_
;
//tmp1 = Gamma(n+alpha+1)/(n! (n+1)^2) dans le cas ou alpha est un entier >0
// = (n+alpha)!/(n! (n+1)^2) c'est le coeff du weight
T
tmp1
=
tN
+
(
T
)
alpha_
;
int
alphaM1
=
alpha_
-
1
;
for
(
int
i
=
1
;
i
<
alphaM1
;
i
++
)
tmp1
*=
(
tmp1
-
(
T
)
1.
);
tmp1
/=
(
tN
+
1.
);
//Approximation by Sould & Secrest
for
(
int
k
=
0
;
k
<
N_
;
k
++
){
if
(
k
==
0
)
{
x
=
(
T
)((
1.0
+
alpha_
)
*
(
3.
+
0.92
*
alpha_
)
/
(
1.
+
2.4
*
N_
+
1.8
*
alpha_
));
}
else
if
(
k
==
1
)
{
x
+=
(
T
)((
15.
+
6.25
*
alpha_
)
/
(
1.
+
0.9
*
alpha_
+
2.5
*
N_
));
}
else
{
//k>=2
int
km1
=
k
-
1
;
T
c0
=
(
T
)((
1.
+
2.55
*
km1
)
/
(
1.9
*
km1
)
+
1.26
*
km1
*
alpha_
/
(
1.
+
3.5
*
km1
));
T
c1
=
(
T
)(
1.
+
0.3
*
alpha_
);
x
+=
c0
*
(
x
-
nodes
(
k
-
2
))
/
c1
;
}
T
xapp
=
x
;
//Refinement
RootUpdate
(
x
);