nan error: star straylight

StrayLight.cpp

@@ -334,6 +334,12 @@ void PointSource(double ju, double sat[3], double ob[3], double py[3], double E[
 
     double E1[7] = { 0, 0, 0, 0, 0, 0, 0 };
     double pstsun, pstmoon;
+    if (isfinite(faisun) == false){
+        faisun = 0;
+    }
+    if (isfinite(thetasun) == false){
+        thetasun = 0;
+    }
     if (thetasun <= PI / 2)
     {
         pstsun = pow(10, PST(faisun, thetasun));
@@ -342,6 +348,12 @@ void PointSource(double ju, double sat[3], double ob[3], double py[3], double E[
     {
         pstsun = 0;
     }
+    if (isfinite(faimoon) == false){
+        faimoon = 0;
+    }
+    if (isfinite(thetamoon) == false){
+        thetamoon = 0;
+    }
     if (thetamoon <= PI / 2 && thetamoon >= PI / 90)
     {
         pstmoon = pow(10, PST(faimoon, thetamoon));
@@ -360,6 +372,7 @@ void PointSource(double ju, double sat[3], double ob[3], double py[3], double E[
         E1[i] = Es_fil[i]*pstsun*Sheltersun*SunRatio[i] + Em_fil[i]*pstmoon*moonphase*Sheltermoon*MoonRatio[i];
     }
     
+    
     for (int i = 0; i < 2; i++)
     {
         ephcom.getPosVel(Planets[i], ju, pos, vel);
@@ -369,6 +382,12 @@ void PointSource(double ju, double sat[3], double ob[3], double py[3], double E[
             continue;
         }
         double faip = OB.Azimuth(planet, PY);
+        if (isfinite(faip) == false){
+            faip = 0;
+        }
+        if (isfinite(thetap) == false){
+            thetap = 0;
+        }
         double pstp = pow(10, PST(faip, thetap));
         double shelterp = ShelterEarth(SAT, planet)*ShelterPlate(planet, OB, PY, fai, op, sum);
         double planetphase = sun.Angle(planet) / PI;
@@ -387,6 +406,12 @@ void PointSource(double ju, double sat[3], double ob[3], double py[3], double E[
             continue;
         }
         double faip = OB.Azimuth(planet, PY);
+        if (isfinite(faip) == false){
+            faip = 0;
+        }
+        if (isfinite(thetap) == false){
+            thetap = 0;
+        }
         double pstp = pow(10, PST(faip, thetap));
         double shelterp = ShelterEarth(SAT, planet)*ShelterPlate(planet, OB, PY, fai, op, sum);
         double proportion = pstp*shelterp/planet.Square();
@@ -418,6 +443,12 @@ void PointSource(double ju, double sat[3], double ob[3], double py[3], double E[
         if (theta<PI / 2 && theta>PI / 90)
         {
             double faistar = OB.Azimuth(star, PY);
+            if (isfinite(faistar) == false){
+                faistar = 0;
+            }
+            if (isfinite(theta) == false){
+                theta = 0;
+            }
             double pst = pow(10, PST(faistar, theta));
             double shelter = ShelterEarth(SAT, star)*ShelterPlate(star,OB,PY,fai,op,sum);//太费时间了
             double proportion = pst*shelter;