diff --git a/ELL1model.C b/ELL1model.C
index 853536c7..6359dd20 100644
--- a/ELL1model.C
+++ b/ELL1model.C
@@ -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)
     {
@@ -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)
diff --git a/readParfile.C b/readParfile.C
index d8a2ae78..ee996ca1 100644
--- a/readParfile.C
+++ b/readParfile.C
@@ -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)
@@ -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);
     }
 }
 
@@ -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 */
@@ -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);
 }
 
 /* ******************************************** */
diff --git a/tests/testTempo2.C b/tests/testTempo2.C
index d12bf7e6..9bf36c72 100644
--- a/tests/testTempo2.C
+++ b/tests/testTempo2.C
@@ -1,6 +1,8 @@
 #include <gtest/gtest.h>
 #include <src/gtest_main.cc>
 #include "tempo2.h"
+#include <string>
+#include <unistd.h>
 #ifndef DATDIR
 #define DATDIR .
 #endif
@@ -8,6 +10,58 @@
 #ifdef LONGDOUBLE_IS_FLOAT128
 #include <quadmath.h>
 #endif
+
+namespace {
+
+const char *ELL1_BASE =
+    "PSRJ J1234+5678\n"
+    "F0 1\n"
+    "PEPOCH 58000\n"
+    "BINARY ELL1\n"
+    "A1 1\n"
+    "TASC 58000\n"
+    "EPS1 0\n"
+    "EPS2 0\n";
+
+const char *DDGR_BASE =
+    "PSRJ J1234+5678\n"
+    "F0 1\n"
+    "PEPOCH 58000\n"
+    "BINARY DDGR\n"
+    "PB 0.5\n"
+    "A1 1\n"
+    "T0 58000\n"
+    "ECC 0.1\n"
+    "OM 45\n"
+    "M2 0.3\n"
+    "MTOT 1.7\n";
+
+void readParText(pulsar *psr, const std::string &text)
+{
+    FILE *fin = tmpfile();
+    ASSERT_TRUE(fin != NULL);
+    ASSERT_NE(fputs(text.c_str(),fin),EOF);
+    rewind(fin);
+    psr->noWarnings = 2;
+    readSimpleParfile(fin,psr);
+    fclose(fin);
+}
+
+double ell1Delay(const std::string &extra, double dt)
+{
+    pulsar psr;
+    MAX_PSR = 1;
+    initialise(&psr,0);
+    readParText(&psr,std::string(ELL1_BASE)+extra);
+    psr.obsn[0].bbat =
+            longdouble(58000.0)+static_cast<longdouble>(dt)/SECDAYl;
+    const double delay = ELL1model(&psr,0,0,-1,0);
+    destroyOne(&psr);
+    return delay;
+}
+
+} // namespace
+
 TEST(testTempo2h, maxpsrset){
     ASSERT_GT(MAX_PSR,0);
 }
@@ -42,6 +96,267 @@ TEST(testLongDouble, printAndParse){
     ASSERT_STREQ("% -1 0.1 2147483649 0.1 2147483649 17179869185 t 50000.1 % x",sb);
 }
 
+TEST(testELL1BinaryPhase, pbFb2MatchesExplicitZeroFb1)
+{
+    const double dt = 1.0e6;
+    const double sparse = ell1Delay(
+            "PB 0.5\n"
+            "FB2 1e-18\n",dt);
+    const double canonical = ell1Delay(
+            "PB 0.5\n"
+            "FB1 0\n"
+            "FB2 1e-18\n",dt);
+    EXPECT_NEAR(sparse,canonical,1e-14);
+}
+
+TEST(testELL1BinaryPhase, sparseFb2ChangesDelay)
+{
+    const double dt = 1.0e6;
+    const double withFb2 = ell1Delay(
+            "PB 0.5\n"
+            "FB2 1e-18\n",dt);
+    const double withoutFb2 = ell1Delay(
+            "PB 0.5\n"
+            "FB2 0\n",dt);
+    EXPECT_GT(fabs(withFb2-withoutFb2),1e-6);
+}
+
+TEST(testELL1BinaryPhase, pbdotSuppliesMissingFb1)
+{
+    const double dt = 1.0e6;
+    const double pbSeconds = 0.5*SECDAY;
+    const double pbdot = 1.0e-8;
+    const double fb1 = -pbdot/(pbSeconds*pbSeconds);
+    // Nested scopes keep only one stack pulsar live at a time:
+    // sizeof(pulsar) is ~5.4 MiB and the default stack is only 8 MiB.
+    const longdouble bbat =
+            longdouble(58000.0)+static_cast<longdouble>(dt)/SECDAYl;
+    double hybridDelay = 0.0;
+    double canonicalDelay = 0.0;
+
+    {
+        pulsar hybrid;
+        MAX_PSR = 1;
+        initialise(&hybrid,0);
+        readParText(&hybrid,std::string(ELL1_BASE)+
+                "PB 0.5\n"
+                "PBDOT 1e-8\n"
+                "FB2 1e-18\n");
+        hybrid.obsn[0].bbat = bbat;
+        hybridDelay = ELL1model(&hybrid,0,0,-1,0);
+        destroyOne(&hybrid);
+    }
+    {
+        pulsar canonical;
+        MAX_PSR = 1;
+        initialise(&canonical,0);
+        readParText(&canonical,std::string(ELL1_BASE)+
+                "PB 0.5\n"
+                "FB1 0\n"
+                "FB2 1e-18\n");
+        canonical.param[param_fb].val[1] = fb1;
+        canonical.obsn[0].bbat = bbat;
+        canonicalDelay = ELL1model(&canonical,0,0,-1,0);
+        destroyOne(&canonical);
+    }
+    EXPECT_NEAR(hybridDelay,canonicalDelay,1e-14);
+}
+
+TEST(testELL1BinaryPhase, fb0PbdotSuppliesMissingFb1)
+{
+    const double dt = 1.0e6;
+    const double fb0 = 2.3148148148148148e-5;
+    const double pbdot = 1.0e-8;
+    const double fb1 = -pbdot*fb0*fb0;
+    const longdouble bbat =
+            longdouble(58000.0)+static_cast<longdouble>(dt)/SECDAYl;
+    double hybridDelay = 0.0;
+    double canonicalDelay = 0.0;
+
+    {
+        pulsar hybrid;
+        MAX_PSR = 1;
+        initialise(&hybrid,0);
+        readParText(&hybrid,std::string(ELL1_BASE)+
+                "FB0 2.3148148148148148e-5\n"
+                "PBDOT 1e-8\n"
+                "FB2 1e-18\n");
+        hybrid.obsn[0].bbat = bbat;
+        hybridDelay = ELL1model(&hybrid,0,0,-1,0);
+        destroyOne(&hybrid);
+    }
+    {
+        pulsar canonical;
+        MAX_PSR = 1;
+        initialise(&canonical,0);
+        readParText(&canonical,std::string(ELL1_BASE)+
+                "FB0 2.3148148148148148e-5\n"
+                "FB1 0\n"
+                "FB2 1e-18\n");
+        canonical.param[param_fb].val[1] = fb1;
+        canonical.obsn[0].bbat = bbat;
+        canonicalDelay = ELL1model(&canonical,0,0,-1,0);
+        destroyOne(&canonical);
+    }
+    EXPECT_NEAR(hybridDelay,canonicalDelay,1e-14);
+}
+
+TEST(testELL1BinaryPhase, missingInteriorFb2IsZero)
+{
+    const double dt = 1.0e5;
+    const double sparse = ell1Delay(
+            "PB 0.5\n"
+            "FB1 -1e-12\n"
+            "FB3 1e-20\n",dt);
+    const double canonical = ell1Delay(
+            "PB 0.5\n"
+            "FB1 -1e-12\n"
+            "FB2 0\n"
+            "FB3 1e-20\n",dt);
+    EXPECT_NEAR(sparse,canonical,1e-14);
+}
+
+TEST(testELL1BinaryPhase, fb0Fb2MatchesEquivalentPbSeries)
+{
+    const double dt = 1.0e6;
+    const double fromFb0 = ell1Delay(
+            "FB0 2.3148148148148148e-5\n"
+            "FB2 1e-18\n",dt);
+    const double fromPb = ell1Delay(
+            "PB 0.5\n"
+            "FB1 0\n"
+            "FB2 1e-18\n",dt);
+    EXPECT_NEAR(fromFb0,fromPb,1e-14);
+}
+
+TEST(testELL1BinaryPhase, ordinaryPbdotPathIsUnchanged)
+{
+    pulsar psr;
+    MAX_PSR = 1;
+    initialise(&psr,0);
+    readParText(&psr,std::string(ELL1_BASE)+
+            "PB 0.5\n"
+            "PBDOT 1e-8\n");
+
+    const double requestedDt = 1.0e6;
+    const double pb = 0.5*SECDAY;
+
+    psr.obsn[0].bbat = longdouble(58000.0)+
+            static_cast<longdouble>(requestedDt)/SECDAYl;
+    const double effectiveDt =
+            (static_cast<double>(psr.obsn[0].bbat)-
+             static_cast<double>(psr.param[param_tasc].val[0]))*SECDAY;
+    const double expectedOrbits =
+            effectiveDt/pb-
+            0.5e-8*(effectiveDt/pb)*(effectiveDt/pb);
+    const double phase = 2.0*M_PI*
+            (expectedOrbits-floor(expectedOrbits));
+    const double dre = sin(phase);
+    const double drep = cos(phase);
+    const double drepp = -sin(phase);
+    const double an = 2.0*M_PI/pb;
+    const double expectedDelay =
+            -dre*(1.0-an*drep+(an*drep)*(an*drep)
+                    +0.5*an*an*dre*drepp);
+
+    EXPECT_NEAR(ELL1model(&psr,0,0,-1,0),expectedDelay,1e-14);
+    destroyOne(&psr);
+}
+
+TEST(testBinaryParfile, explicitFb1DisablesPbdot)
+{
+    pulsar psr;
+    MAX_PSR = 1;
+    initialise(&psr,0);
+    readParText(&psr,std::string(ELL1_BASE)+
+            "PB 0.5\n"
+            "PBDOT 1e-12 1\n"
+            "FB1 -1e-18\n");
+
+    EXPECT_EQ(psr.param[param_fb].paramSet[1],1);
+    EXPECT_EQ(psr.param[param_pbdot].paramSet[0],0);
+    EXPECT_EQ(psr.param[param_pbdot].fitFlag[0],0);
+    destroyOne(&psr);
+}
+
+TEST(testBinaryParfile, pbdotRemainsActiveWithoutExplicitFb1)
+{
+    pulsar psr;
+    MAX_PSR = 1;
+    initialise(&psr,0);
+    readParText(&psr,std::string(ELL1_BASE)+
+            "PB 0.5\n"
+            "PBDOT 1e-12 1\n"
+            "FB2 1e-18\n");
+
+    EXPECT_EQ(psr.param[param_fb].paramSet[1],0);
+    EXPECT_EQ(psr.param[param_pbdot].paramSet[0],1);
+    EXPECT_EQ(psr.param[param_pbdot].fitFlag[0],1);
+    destroyOne(&psr);
+}
+
+TEST(testELL1BinaryPhase, explicitFb1HasZeroPbdotDerivative)
+{
+    pulsar psr;
+    MAX_PSR = 1;
+    initialise(&psr,0);
+    readParText(&psr,std::string(ELL1_BASE)+"PB 0.5\n");
+
+    // Bypass post-parse consistency handling to exercise the defensive
+    // derivative branch for direct programmatic construction.
+    psr.param[param_pbdot].paramSet[0] = 1;
+    psr.param[param_pbdot].val[0] = 1e-12;
+    psr.param[param_fb].paramSet[1] = 1;
+    psr.param[param_fb].val[1] = -1e-18;
+    psr.obsn[0].bbat =
+            longdouble(58000.0)+longdouble(1.0e5)/SECDAYl;
+
+    EXPECT_EQ(ELL1model(&psr,0,0,param_pbdot,0),0.0);
+    destroyOne(&psr);
+}
+
+TEST(testBinaryParfile, ddgrRejectsFb0)
+{
+    ::testing::FLAGS_gtest_death_test_style = "threadsafe";
+    EXPECT_EXIT(
+        {
+            ASSERT_EQ(dup2(STDERR_FILENO,STDOUT_FILENO),STDOUT_FILENO);
+            pulsar psr;
+            MAX_PSR = 1;
+            initialise(&psr,0);
+            readParText(&psr,std::string(DDGR_BASE)+"FB0 2e-5\n");
+        },
+        ::testing::ExitedWithCode(1),
+        "BINARY DDGR does not support FB parameters");
+}
+
+TEST(testBinaryParfile, ddgrRejectsHigherFb)
+{
+    ::testing::FLAGS_gtest_death_test_style = "threadsafe";
+    EXPECT_EXIT(
+        {
+            ASSERT_EQ(dup2(STDERR_FILENO,STDOUT_FILENO),STDOUT_FILENO);
+            pulsar psr;
+            MAX_PSR = 1;
+            initialise(&psr,0);
+            readParText(&psr,std::string(DDGR_BASE)+"FB1 -1e-18\n");
+        },
+        ::testing::ExitedWithCode(1),
+        "BINARY DDGR does not support FB parameters");
+}
+
+TEST(testBinaryParfile, ddgrWithoutFbParses)
+{
+    pulsar psr;
+    MAX_PSR = 1;
+    initialise(&psr,0);
+    readParText(&psr,DDGR_BASE);
+    EXPECT_STREQ(psr.binaryModel,"DDGR");
+    EXPECT_EQ(psr.param[param_fb].paramSet[0],0);
+    EXPECT_EQ(psr.param[param_fb].paramSet[1],0);
+    destroyOne(&psr);
+}
+
 TEST(testFormResiduals, basicBATs){
     pulsar _psr;
     pulsar *psr = &_psr;