Normalization for projection matrices?

[Leo Stein]

Hi,

I'm interested in calculating multipole moments from stellar models computed by ESTER. For example, I'd like to compute the integrals

  M_\ell = 2 \pi \iint \rho(r,\theta) P_\ell (\mu) d\mu r^{\ell +2} dr

for even \ell, where \mu = \cos\theta, and another set of integrals

  S_\ell = \frac{4 \pi}{\ell+1}  \iint \rho(r,\theta)  \Omega(r,\theta) \frac{d P_\ell (\mu)}{d\mu} \sin^2\theta  d\mu r^{\ell +3} dr

for odd \ell (image of these equations attached).

Now these integrals already involve performing a decomposition in terms of Legendre polynomials, so I thought this should be quite straightforward with ESTER. This is what I tried for the M_l integrals:

star2d s;

s.read(argv[1]);

matrix rho_coeffs = (s.rho, s.map.leg.P_00);

matrix ell_vect = vector(0.0, (double)s.nth - 1.0, s.nth);

matrix ells(s.nr, s.nth);

for ( int i=0; i < s.nr; i++)

  ells.setrow(i, ell_vect);

// We want weights of r^(\ell + 2)

matrix rweights = pow( s.r, 2.0 + ells );

// Don't forget the Jacobian factor coming from dr/dzeta

matrix multipoles = 2*PI*(s.map.gl.I, rho_coeffs * rweights * s.map.rz);

matrix unit_conv = s.units.rho * pow(s.units.r, 3.0 + ell_vect );

multipoles *= unit_conv; 

multipoles.write();

Here I am trying to use P_00 to convert the data (for rho) from the collocation representation to the harmonic coefficients rho_ell. Am I correct in assuming that column number \ell of rho_coeffs now contains \int \rho P_\ell(\mu) d\mu ? Or is there a normalization factor?

The zeroth moment should correspond to the total mass of the star. However, the above code gives a number which is different by a factor of 1.98. This is why I suspect a normalization issue.

I haven't yet tried to compute the S_l integrals because I am not sure that I know the right basis functions. I am interested in odd \ell and finding coefficients of dP_\ell, so is it correct to use P_10 to extract the harmonic coefficients?

Thanks in advance,

Leo

[Paco Espinosa]

Hi Leo,

Yes, there is a normalization factor. The normalized Legendre polynomials are defined in ESTER as

Pnorm_l = sqrt(2*l+1) * P_l

such that the scalar product between two norm. Leg. polynomials is 2*delta(l,l').

However, I've found a problem in your code, given that both r and dr/dzeta also depend on theta. Personally, I will choose a more direct method. In
ESTER, you have the following matrices:

map.leg.P1_00 : Row n is Pnorm_(2n)  (starting with n=0)
map.leg.P1_01 : Row n is Pnorm_(2n+1)
map.leg.P1_10 : Row n is dPnorm_(2n+1)/dtheta
map.leg.P1_11 : Row n is dPnorm_(2n+2)/dtheta

Then, to calculate Ml, I would do something like:

matrix Pl = s.map.leg.P1_00.row(l/2) / sqrt(2*l+1);
double Ml = 2*PI * (s.map.gl.I, s.rho * Pl * pow(s.r,l+2) * s.map.rz, s.map.leg.I_00)(0);
Ml *= s.units.rho * pow(s.units.r, l+3);

To finish, just a little advice. The models are not spherically symmetric, that means that the decomposition in spherical harmonics is not unique.
It depends on the mapping and more precisely on the way that we choose the surfaces zeta=constant.  So the code above may lead to
slightly different results if you use a different number of domains or a different radial resolution. I let you think about that.

Best,

Paco

--
You received this message because you are subscribed to the Google Groups "ester-project" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ester-projec...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Paco Espinosa]

Hi Leo,

just forget my last comment. It is true that the spectral decomposition won't be the same for different mappings, but the result of the integral (Ml) should
be the same independently of how the mapping is chosen.

Paco.

[Michel Rieutord]

Le 16/04/2014 16:40, Paco Espinosa a écrit :

Hi Leo,

just forget my last comment. It is true that the spectral decomposition won't be the same for different mappings, but the result of the integral (Ml) should
be the same independently of how the mapping is chosen.

I agree. I was surprised of your remark, Paco (I thought I overlooked something!).

Michel

[Leo Stein]

Dear Paco,

Thanks for clarifying the normalization conventions used in the code. If you understood what I was trying to do in my original code, you saw that I was trying to compute all multipole moments at the same time instead of one by one (thinking this would be more efficient). This is actually not a very expensive calculations so I instead integrate (as you recommended) to get one multipole moment at a time.

In case this is useful for anybody in the future, I'm including these utility functions for calculating mass and mass-current moments from a stellar model in code units, cgs, and geometric units. By the way, can you clarify what the internal code units are for r, rho, omega?

Best

Leo

---

double M(star2d &s, unsigned int ell)

{

  if (1 == ell%2)

    return 0;

  matrix Pl = s.map.leg.P1_00.row(ell/2) / sqrt(2*ell+1);

  double Ml = 2.0*PI * (s.map.gl.I,

                      s.rho * Pl * pow(s.r,ell+2.0) * s.map.rz,

                      s.map.leg.I_00)(0);

  return Ml;

};

double Mcgs(star2d &s, unsigned int ell)

{

  double Ml = M(s, ell);

  Ml *= s.units.rho * pow(s.units.r, ell+3.0);

  return Ml;

};

double Mgeom(star2d &s, unsigned int ell)

{

  double Ml = Mcgs(s, ell);

  Ml *= GRAV / (C_LIGHT * C_LIGHT);

  return Ml;

};

double S(star2d &s, unsigned int ell)

{

  if (0 == ell%2)

    return 0;

  matrix dPl = s.map.leg.P1_10.row((ell-1)/2) / sqrt(2.0*ell+1.0);

  double Sl = - 4.0*PI / (ell + 1.0) *

               (s.map.gl.I,

                s.rho * s.w * sin(s.th)

                * dPl * pow(s.r,ell+3.0) * s.map.rz,

                s.map.leg.I_00)(0);

  return Sl;

};

double Scgs(star2d &s, unsigned int ell)

{

  double Sl = S(s, ell);

  Sl *= s.units.rho * s.units.Omega * pow(s.units.r, ell+4.0);

  return Sl;

};

double Sgeom(star2d &s, unsigned int ell)

{

  double Sl = Scgs(s, ell);

  Sl *= GRAV / pow(C_LIGHT, 3.0);

  return Sl;

};

[Paco Espinosa]

The units used internally by ESTER are normalized using the polar radius (R), central  density (rhoc), central pressure (pc) and central temperature (Tc).

Units:
- Density: rhoc
- Pressure: pc
- Temperature: Tc
- Length: R
- Omega: sqrt(pc/rhoc)/R

Note that this only concerns the internal calculations in the code and could change in future versions. The safest way to work with this units is to
use the units structure defined inside the class star2d (and star1d).

class star2d {

[...]

struct units_struct {
   double rho,p,phi,T,Omega,r,v,F;
} units;

[...]

};

where rho->density, p->pressure, phi->(gravitational) potential, T->temperature, Omega->Angular velocity, r->length, v->velocity
and F->acceleration (specific force). These values are used to convert code units to cgs.

Best,
Paco