Added condition on xa>xb in tides calculation

src/tides/tidal_lambda2.c

@@ -156,34 +156,54 @@ double tidal_lambda2(struct stardata_t * const stardata,
                                                  RHS);
                 Dprint("Check func at root xa = %g : %g\n",
                        xa,Zahn_Eq_14b(xa,RHS));
+                Dprint("Check xa falls into CZ with xa>xb xa = %g, xb %g\n",
+                       xa,xb);
 
-                const double integral1 = GSL_integrator(xa, // lower
-                                                        1.0, // upper
-                                                        integration_tolerance,
-                                                        &integration_error,
-                                                        GSL_INTEGRATOR_QAG,
-                                                        &Zahn_integral_13);
+                if(xa>xb)
+                {
+                    const double integral1 = GSL_integrator(xa, // lower
+                                                            1.0, // upper
+                                                            integration_tolerance,
+                                                            &integration_error,
+                                                            GSL_INTEGRATOR_QAG,
+                                                            &Zahn_integral_13);
 
-                const double integral2 = GSL_integrator(xb, // lower
-                                                        xa, // upper
-                                                        integration_tolerance,
-                                                        &integration_error,
-                                                        GSL_INTEGRATOR_QAG,
-                                                        &Zahn_integral_14a)
-                    * pow(xa,7.0/6.0) * Pow1p5(1.0 - xa);
+                    const double integral2 = GSL_integrator(xb, // lower
+                                                            xa, // upper
+                                                            integration_tolerance,
+                                                            &integration_error,
+                                                            GSL_INTEGRATOR_QAG,
+                                                            &Zahn_integral_14a)
+                       * pow(xa,7.0/6.0) * Pow1p5(1.0 - xa);
 
-                integral = integral1 + integral2;
-                Dprint("have roots: xb = %g, xa = %g, integral = %g + %g = %g\n",
+                    integral = integral1 + integral2;
+                    Dprint("have roots: xb = %g, xa = %g, integral = %g + %g = %g\n",
+                           xb,
+                           xa,
+                           integral1,
+                           integral2,
+                           integral);
+                }
+                else
+                {
+                /*
+                 * There is a root but if falls below convective zone,  single integral in Eq. 13
+                 */
+                integral = GSL_integrator(xb, // lower
+                                          1.0, // upper
+                                          integration_tolerance,
+                                          &integration_error,
+                                          GSL_INTEGRATOR_QAG,
+                                          &Zahn_integral_13);
+                Dprint("irrelevant root ignored: xb = %g, integral = %g\n",
                        xb,
-                       xa,
-                       integral1,
-                       integral2,
                        integral);
+                }
             }
             else
             {
                 /*
-                 * No roots, single integral in Eq. 13
+                 * No roots,  single integral in Eq. 13
                  */
                 integral = GSL_integrator(xb, // lower
                                           1.0, // upper