inf (divide-by-zero) in NS::CalculateAndSetPulsarParameters, NS::CalculateSpinDownRate(), and NS::SpinDownIsolatedPulsar()

[jeffriley]

Running COMPAS with:

./compas -n1 --random-seed 0

results in inf in NS::CalculateAndSetPulsarParameters in the following statement:

m_PulsarDetails.spinFrequency = _2_PI / (m_PulsarDetails.spinPeriod * SECONDS_IN_MS);

because CalculateBirthSpinPeriod() in this statement:

m_PulsarDetails.spinPeriod = CalculateBirthSpinPeriod();                // spin period in ms

returns 0.0 because the default PULSAR_BIRTH_SPIN_PERIOD_DISTRIBUTION is ZERO

void NS::CalculateAndSetPulsarParameters() {

    m_PulsarDetails.magneticField     = PPOW(10.0, CalculateBirthMagneticField()) * GAUSS_TO_TESLA;         // magnetic field in Gauss -> convert to Tesla
    m_PulsarDetails.spinPeriod        = CalculateBirthSpinPeriod();                                         // spin period in ms

    m_PulsarDetails.spinFrequency     = _2_PI / (m_PulsarDetails.spinPeriod * SECONDS_IN_MS);
    m_PulsarDetails.birthPeriod       = m_PulsarDetails.spinPeriod * SECONDS_IN_MS;                         // convert from ms to s 
    
    m_MomentOfInertia_CGS             = CalculateMomentOfInertiaCGS();                                      // in CGS g cm^2
	
    // Note we convert neutronStarMomentOfInertia from CGS to SI here
    m_PulsarDetails.spinDownRate      = CalculateSpinDownRate(m_PulsarDetails.spinFrequency, m_MomentOfInertia_CGS, 
    m_PulsarDetails.magneticField, m_Radius * RSOL_TO_KM);  
    m_PulsarDetails.birthSpinDownRate = m_PulsarDetails.spinDownRate; 
    m_AngularMomentum_CGS             = m_MomentOfInertia_CGS * m_PulsarDetails.spinFrequency;              // in CGS g cm^2 s^-1
}

The naive solution is to change the statement:

m_PulsarDetails.spinFrequency = _2_PI / (m_PulsarDetails.spinPeriod * SECONDS_IN_MS);

to:

m_PulsarDetails.spinFrequency = m_PulsarDetails.spinPeriod > 0.0 ? _2_PI / (m_PulsarDetails.spinPeriod * SECONDS_IN_MS) : 0.0;   // don't use utils::Compare() here

so that m_PulsarDetails.spinFrequency is set to 0.0 if CalculateBirthSpinPeriod() returns 0.0

@yuzhesong , @SimonStevenson, can you confirm that the suggestion above is the best solution?

A related problem occurs in NS::CalculateSpinDownRate(), when parameter p_Omega is passed as 0.0. This statement:

double period = _2_PI / p_Omega;                        // convert frequency to period

results in inf

double NS::CalculateSpinDownRate(const double p_Omega, const double p_MomentOfInteria, const double p_MagField, const double p_Radius) const {

    // pow() is slow - use multiplication

    double period            = _2_PI / p_Omega;                                                                 // convert frequency to period
    double cgsRadius         = p_Radius * KM_TO_CM;                                                             // radius in cm
    double radius_6          = cgsRadius * cgsRadius * cgsRadius * cgsRadius * cgsRadius * cgsRadius;
    double cgsMagField       = p_MagField * TESLA_TO_GAUSS;                                                     // B field in G
    double magField_2        = cgsMagField * cgsMagField;
    constexpr double _8_PI_2 = 8.0 * PI_2;
    constexpr double _3_C_3  = 3.0E6 * C * C * C;                                                               // 3.0 * (C * 100.0) * (C * 100.0) * (C * 100.0)
    double pDotTop           = _8_PI_2 * radius_6 * magField_2;
    double pDotBottom        = _3_C_3 * p_MomentOfInteria * period;
    double pDot              = pDotTop / pDotBottom;                                                            // period derivative 
   
    return -pDot * p_Omega / period;                                                                            // convert period derivative to frequency derivative, which is what is recorded in the output
}

The naive solution is to add the following sanity check as the first line of the function:

if (p_Omega<= 0.0) return 0.0;                                  // don't use utils::Compare() here

so we just return a spin down rate of 0.0 if the spin freqency passed in p_Omega is 0.0

@yuzhesong , @SimonStevenson, can you confirm that the suggestion above is the best solution?

A further related problem becomes evident when the above two problems are addressed: the following statement in NS::SpinDownIsolatedPulsar() results in inf:

double initialSpinPeriod = _2_PI / m_PulsarDetails.spinFrequency;

which can naively be fixed by replacing the statement with:

double initialSpinPeriod = m_PulsarDetails.spinFrequency > 0.0 ? _2_PI / m_PulsarDetails.spinFrequency : 0.0;    // don't use utils::Compare() here

But maybe there's a better solution - can we avoid calling some of these functions if the spin frequency is 0.0?

@yuzhesong, @SimonStevenson Is there a better fix overall when spin frequency is 0.0?

To Reproduce
Run COMPAS with:

./compas -n1 --random-seed 0

Expected behavior
Production code does not cause inf

**Versioning

OS: Ubuntu 22.04
COMPAS v03.07.01

[SimonStevenson]

I mean, the frequency really is infinite when the birth spin period is zero (non rotating). It's not a physically very sensible choice, but it is correct. I don't think setting them to 0 is the right thing to do (that corresponds to an infinite spin period, which is not consistent).

The spin down rate should just return 0 if the spin period is 0 (it can't spin down any more than non-spinning).

I wouldn't really call this a major-bug either.

[ilyamandel]

As an aside, why do we need to store m_PulsarDetails.spinFrequency in memory if we can trivially compute it as a division?

[jeffriley]

A couple of things:

the frequency really is infinite when the birth spin period is zero (non rotating)

How can the frequency with which a body is rotating be infinite if it isn't rotating? Surely if something isn't rotating, the frequency with which it is rotating is zero. If anything should be infinite shouldn't it be the period (the time it takes for one full rotation)?

I wouldn't really call this a major-bug either.

The problem is that we divide-by-zero. Not only shouldn't we do that, if I turn fp-error handling on (not available yet on macOS), I don't get past this divide-by-zero... so fp-error handling is rendered almost useless...

[ilyamandel]

Right, spin period = 0 is a bit of a special case. As spin periods get lower and lower, the frequency gets higher and higher.
Technically, a spin period of zero, if taken as a limit, should correspond to infinitely fast rotation. Thinking further (and correcting a comment I made above), we should in fact store spin frequencies, not spin periods, and work with these only. A spin frequency of zero is unambiguously a sign of no rotation (as frequencies decrease, the body spins slower, so this is a natural limit of infinitely slow rotation).

[jeffriley]

Right, but here it's not the limit - a birth spin period of 0.0 here is meant to indicate it's not spinning (as you noted above ('...when the birth spin period is zero (non rotating)')).

we should in fact store spin frequencies, not spin periods, and work with these only.
A spin frequency of zero is unambiguously a sign of no rotation

I agree with that - let's do that :-) And if we need to calculate the period from frequency we should nt divide by zero - if the star isn't spinning we could avoid a lot of the calculation in the NS code, couldn't we?

[ilyamandel]

(as you noted above ('...when the birth spin period is zero (non rotating)')).

That was @SimonStevenson , not me [just giving credit :) , I agree with your solution]

[ilyamandel]

@yuzhesong, I think we agreed to switch to storing spin frequencies, rather than spin periods, to avoid the ambiguity of "zero period" and the associated issues. Can you please implement this as part of the PR you are currently working on?