@@ -0,0 +1,2116 @@
+# coding: utf-8%0A%0A%22%22%22 Check how timestep and number of steps for an orbit effects the frequencies. %22%22%22%0A%0Afrom __future__ import division, print_function%0A%0A__author__ = %22adrn %3Cadrn@astro.columbia.edu%3E%22%0A%0A# Standard library%0Aimport os%0Aimport cPickle as pickle%0A%0A# Third-party%0Aimport numpy as np%0A%0A# Project%0Aimport gary.dynamics as gd%0Aimport gary.integrate as gi%0A%0A# project%0Afrom streammorphology.potential import potential_registry%0A%0Adef main():%0A    filename = %22/home/adrian/freq_dicts.pickle%22%0A    potential = potential_registry%5B'triaxial-NFW'%5D%0A%0A    # slow and fast diffusing orbits%0A    w0 = np.array(%5B%5B-17.504234723,-17.2283745157,-9.07711397761,0.0721992194851,0.021293129758,0.199775306493%5D,%0A                   %5B-1.69295332221,-13.78418595,15.6309115075,0.0580704842,0.228735516722,0.0307028904261%5D%5D)%0A%0A    dts = np.array(%5B0.1, 0.25, 0.5, 1.%5D)%0A    t2 = 100000%0A    nsteps = t2 / dts%0A    dts,nsteps = map(np.ravel, np.meshgrid(dts, nsteps))%0A%0A    if not os.path.exists(filename):%0A        freq_dicts = %5Bdict(), dict()%5D%0A        for dt,nstep in zip(dts,nsteps):%0A            print(dt,nstep)%0A            t,w = potential.integrate_orbit(w0, dt=dt, nsteps=nstep,%0A                                            Integrator=gi.DOPRI853Integrator)%0A%0A            for i in range(2):%0A                print(%22%5Ct%7B%7D%22.format(i))%0A                naff = gd.NAFF(t)%0A                fs = gd.poincare_polar(w%5B:,i%5D)%0A                try:%0A                    freqs,d,nvec = naff.find_fundamental_frequencies(fs)%0A                except:%0A                    print(%22Orbit %7B%7D failed for dt=%7B%7D, nsteps=%7B%7D%22.format(i,dt,nstep))%0A                    continue%0A                freq_dicts%5Bi%5D%5B(dt,nstep)%5D = list(freqs)%0A%0A        with open(filename, %22w%22) as f:%0A            pickle.dump(freq_dicts, f)%0A%0A    with open(filename, %22r%22) as f:%0A        freq_dicts = pickle.load(f)%0A%0A    for i in range(2):%0A        print(%22-----------------------------------%22)%0A        print(%22Orbit %7B%7D%22.format(i))%0A        base_freqs = np.array(freq_dicts%5Bi%5D%5B(0.1,1000000)%5D)%0A        for k,v in freq_dicts%5Bi%5D.items():%0A            print(k, np.array(v) - base_freqs)%0A%0Aif __name__ == '__main__':%0A    main()%0A%0A

add script