Skip to content

Evaluate sparse ELL1 FB derivatives and reject inactive DDGR FB terms #7

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
vhaasteren wants to merge 1 commit into from
Closed

Evaluate sparse ELL1 FB derivatives and reject inactive DDGR FB terms #7

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
50 changes: 36 additions & 14 deletions ELL1model.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -112,24 +112,42 @@ double ELL1model(pulsar *psr,int p,int ipos,int param,int k)

ct = psr[p].obsn[ipos].bbat;
tt0 = (ct-t0asc)*SECDAY;
// --- Changes to handle higher orbital-frequency derivatives (FB1, FB2, ...) ---
// 04/2015, H. J. Pletsch
// --- Handle higher orbital-frequency derivatives (FB1, FB2, ...) ---
// A missing coefficient is zero; it does not disable higher terms.
int useHigherFB = 0;
for (int j=1;j<psr[p].param[param_fb].aSize;j++)
{
if (psr[p].param[param_fb].paramSet[j]==1)
{
useHigherFB = 1;
break;
}
}

orbits = tt0/pb;
if (psr[p].param[param_fb].paramSet[1]==1) {
int j;
if (useHigherFB == 1)
{
double fac = 1.0;
for (j=1;j<psr[p].param[param_fb].aSize;j++) {
double fbx;
fac = fac/((double)(j+1));
if (psr[p].param[param_fb].paramSet[j]==1) {
fbx = psr[p].param[param_fb].val[j];
orbits += fac * fbx * pow(tt0,j+1);
}
for (int j=1;j<psr[p].param[param_fb].aSize;j++)
{
fac = fac/((double)(j+1));
if (psr[p].param[param_fb].paramSet[j]==1)
{
const double fbx = psr[p].param[param_fb].val[j];
orbits += fac*fbx*pow(tt0,j+1);
}
else if (j==1 && psr[p].param[param_pbdot].paramSet[0]==1)
{
const double fb1 = -pbdot/(pb*pb);
orbits += fac*fb1*pow(tt0,2);
}
}
} else {
}
else
{
orbits -= 0.5*(pbdot+xpbdot)*pow(tt0/pb,2);
}
// --- End of changes to handle higher orbital-frequency derivatives ---
}
// --- End higher orbital-frequency derivatives ---

if (psr[p].param[param_orbifunc].paramSet[0] == 1)
{
Expand Down Expand Up @@ -240,7 +258,11 @@ double ELL1model(pulsar *psr,int p,int ipos,int param,int k)
else if (param==param_eps2dot)
return Ceps2*tt0;
else if (param==param_pbdot)
{
if (psr[p].param[param_fb].paramSet[1]==1)
return 0.0;
return 0.5*tt0*(-Csigma*an*SECDAY*tt0/(pb*SECDAY));
}
else if (param==param_a1dot)
return Cx*tt0;
else if (param==param_sini)
Expand Down
35 changes: 35 additions & 0 deletions readParfile.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,7 @@ void getValue(char *str,int v1,int v2,pulsar *psr,int l,int arr);
void removeCR(char *str);
void checkLine(pulsar *p,char *str,FILE *fin,parameter *elong,parameter *elat);
void checkAllSet(pulsar *psr,parameter elong,parameter elat,char *filename);
static void checkBinaryParameterConsistency(pulsar *psr);

/* Function to set up default parameters before reading a .par file */
int setupParameterFileDefaults(pulsar *psr)
Expand Down Expand Up @@ -150,6 +151,7 @@ void readParfileGlobal(pulsar *psr,int npsr,char tpar[MAX_STRLEN][MAX_FILELEN],
checkLine(psr+p,str,fin,&elong,&elat);
}
fclose(fin);
checkBinaryParameterConsistency(psr+p);
}
}

Expand Down Expand Up @@ -2269,6 +2271,38 @@ void checkLine(pulsar *psr,char *str,FILE *fin,parameter *elong, parameter *elat
}
}

static void checkBinaryParameterConsistency(pulsar *psr)
{
int anyFB = 0;
for (int j=0;j<psr->param[param_fb].aSize;j++)
{
if (psr->param[param_fb].paramSet[j]==1)
{
anyFB = 1;
break;
}
}

if (strcasecmp(psr->binaryModel,"DDGR")==0 && anyFB==1)
{
logerr("BINARY DDGR does not support FB parameters; "
"use PB/PBDOT or select a binary model that evaluates FBn.");
exit(1);
}

if (strcasecmp(psr->binaryModel,"ELL1")==0 &&
psr->param[param_fb].paramSet[1]==1 &&
psr->param[param_pbdot].paramSet[0]==1)
{
displayMsg(1,"BINFB1",
"Both PBDOT and FB1 are set for ELL1; explicit FB1 takes "
"precedence and PBDOT is disabled.",
"",psr->noWarnings);
psr->param[param_pbdot].paramSet[0] = 0;
psr->param[param_pbdot].fitFlag[0] = 0;
}
}

void checkAllSet(pulsar *psr,parameter elong,parameter elat,char *filename)
{
/* Check if we have read a position epoch */
Expand Down Expand Up @@ -2384,6 +2418,7 @@ void checkAllSet(pulsar *psr,parameter elong,parameter elat,char *filename)
printf("-> Using %s instead.(You should set CLK to what you really mean!)\n", clk);
strcpy(psr->clock, clk);
}
checkBinaryParameterConsistency(psr);
}

/* ******************************************** */
Expand Down
Loading