Eccentric anomaly failed to converge

[astropatel]

Dear all,

I encountered this "failure to converge" error while trying to run :

res = sim.SurveySimulation.run_sim()

using the "template_WFIRST_EarthTwinHabZone.json" scriptfile.

It seems to have failed after the last mission times listed below:

Current mission time:  0.0 d
Current mission time:  9.60003167231 d
Current mission time:  11.6313483646 d
Current mission time:  91.0376767603 d
Current mission time:  92.0376767603 d
Current mission time:  93.0376767603 d
Current mission time:  173.99003259 d
Current mission time:  175.073673146 d
Current mission time:  262.617593242 d
Current mission time:  342.664090476 d
Current mission time:  344.424390472 d
Current mission time:  345.424390472 d
Current mission time:  346.424390472 d
Current mission time:  347.424390472 d
Current mission time:  424.387741558 d

I looked at the while loop in eccanom.py where the convergence occurs and printed out the errors the piece calculates as part of the condition for convergence. It seems the eccentric anomaly converges after only a few iterations as shown below. The entire run of 200 iterations produces 1.1368e-13 by the end of the run:

err_arr = [1.0, 48.815465289295275, 8.0299584348608732, 2.0708116579178295, 0.085908077132557992, 0.00048370889385296323, 1.4405713955056854e-08, 1.1368683772161603e-13, 1.1368683772161603e-13, ..., 
1.1368683772161603e-13]

However, it doesn't ever reduce to the required tolerance level set by

tolerance = np.finfo(float).eps*4.01

which, in my case = 8.903988657493755e-16.

Any ideas? I've listed the full traceback below.

Exception                                 Traceback (most recent call last)
<ipython-input-112-667210cb35ae> in <module>()
----> 1 res = sim.SurveySimulation.run_sim()

/Users/rpatel/Dropbox/Research/WFIRST/EXOSIMS/EXOSIMS/Prototypes/SurveySimulation.pyc in run_sim(self)
    300 
    301             # acquire a new target star index
--> 302             sInd, DRM = self.next_target(sInd, DRM)
    303 
    304             # append result values to self.DRM

/Users/rpatel/Dropbox/Research/WFIRST/EXOSIMS/EXOSIMS/Prototypes/SurveySimulation.pyc in next_target(self, sInd, DRM)
    342         while not TK.mission_is_over():
    343             # find keepout Boolean values (kogood)
--> 344             Obs.keepout(TK.currentTimeAbs, TL, OS.telescopeKeepout)
    345             # if observable targets, pick one, else allocate time and try again
    346             if np.any(Obs.kogood):

/Users/rpatel/Dropbox/Research/WFIRST/EXOSIMS/EXOSIMS/Observatory/WFIRSTObservatory.pyc in keepout(self, time, catalog, koangle)
     92                     self.solarSystem_body_position(time, 'Venus'), # Venus
     93                     self.solarSystem_body_position(time, 'Earth'), # Earth
---> 94                     self.solarSystem_body_position(time, 'Mars'), # Mars
     95                     self.solarSystem_body_position(time, 'Jupiter'), # Jupiter
     96                     self.solarSystem_body_position(time, 'Saturn'), # Saturn

/Users/rpatel/Dropbox/Research/WFIRST/EXOSIMS/EXOSIMS/Prototypes/Observatory.pyc in solarSystem_body_position(self, time, bodyname)
    307             return self.spk_body(time,bodyname)
    308         else:
--> 309             return self.keplerplanet(time,bodyname)
    310 
    311     def spk_body(self, time, bodyname):

/Users/rpatel/Dropbox/Research/WFIRST/EXOSIMS/EXOSIMS/Prototypes/Observatory.pyc in keplerplanet(self, time, bodyname)
    408         wp = w - O
    409         # Find eccentric anomaly
--> 410         E = eccanom(M,e)[0]
    411         # Find true anomaly
    412         nu = np.arctan2(np.sin(E) * np.sqrt(1 - e**2), np.cos(E) - e)

/Users/rpatel/Dropbox/Research/WFIRST/EXOSIMS/EXOSIMS/util/eccanom.pyc in eccanom(M, e)
     54 
     55     if numIter == maxIter:
---> 56         raise Exception("eccanom failed to converge.")
     57 
     58     return E

Exception: eccanom failed to converge.

[dsavransky]

I'll take a look. It might be that we just need to relax this tolerance. Note, however, that you are going down this code branch because you do not have jplephem installed. Using that module should produce more robust (and faster) performance for solar system bodies.

[dsavransky]

I finally identified the root cause of this failure - the ephemeris propagation was generating mean anomaly values outside of [0,2*pi), which caused the Newton-Raphson to freak out. Fixed by adding appropriate mod statements to the Observatory prototype, and verified that eccanom never fails to converge for any value of M in the valid range.