#!/usr/bin/env python

from math import *

R0 = 8500.
Theta0 = 220.
tend = 1*2*pi*R0/Theta0
dt = 0.0001
dtp = 1.
R10 = 500.


def F(x, y):
    R = sqrt(x*x+y*y)
    phi = atan2(y, x)
    Fr = Theta0**2/R
    return (-cos(phi)*Fr, -sin(phi)*Fr)


x = R0+R10
y = 0
vx = 0
h0 = Theta0*R0
vy = -h0/(R0+R10)
kappa = sqrt(3.*h0**2/R0**4 - Theta0**2/R0**2)
#print h0, R0, Theta0, kappa, kappa*R0
t = 0
tp = 0
ip = 0

fgp = file('star2.gp', 'w')
fgp.write('set terminal postscript eps enhanced solid color 15\n')
fgp.write('set parametric\n')
fgp.write('set size ratio -1\n')
fgp.write('set label 10 at 0,0 point ps 1\n')

print '% 10.3e % 10.3e % 10.3e % 10.3e % 10.3e' \
% (t, x, y, vx, vy)
while t < tend:

    if (t-tp) > dtp:
        ip += 1
        tp = t
        print '% 10.3e % 10.3e % 10.3e % 10.3e % 10.3e' \
        % (t, x, y, vx, vy)
        fgp.write('set output "images2/star%04d.eps"\n' % (ip))        
        fgp.write('set label 11 at %12.5e, %12.5e point ps 2 pt 3\n' \
        % (x, y))        
        fgp.write('set label 12 at %12.5e, %12.5e point ps 1 pt 2\n' \
        % (R0*cos(-Theta0*t/R0), R0*sin(-Theta0*t/R0)))        
        phi = atan2(y, x)
        fgp.write(('plot "star.dat" u 2:3 ev ::::%d t "" w l,' 
        + ' %12.5e*sin(t), %12.5e*cos(t) t "", '\
        + '%12.5e*sin(t)+%12.5e*cos(t) + %12.5e'\
        ', %12.5e*sin(t)+%12.5e*cos(t) + %12.5e t ""\n') \
        % (ip, R0, R0 \
        , R10*cos(phi), -2*Theta0*R10/R0/kappa*sin(phi), R0*cos(-Theta0*t/R0) \
        , R10*sin(phi), 2*Theta0*R10/R0/kappa*cos(phi), R0*sin(-Theta0*t/R0)))        
        fgp.write('unset label 11\n')        
        fgp.write('unset label 12\n')        
    
    dvx1 = F(x, y)[0]*dt
    dvy1 = F(x, y)[1]*dt
    dx1 = vx*dt
    dy1 = vy*dt

    dvx2 = F(x+0.5*dx1, y+0.5*dx1)[0]*dt
    dvy2 = F(x+0.5*dx1, y+0.5*dx1)[1]*dt
    dx2 = (vx+0.5*dvx1)*dt
    dy2 = (vy+0.5*dvy1)*dt

    dvx3 = F(x+0.5*dx2, y+0.5*dx2)[0]*dt
    dvy3 = F(x+0.5*dx2, y+0.5*dx2)[1]*dt
    dx3 = (vx+0.5*dvx2)*dt
    dy3 = (vy+0.5*dvy2)*dt

    dvx4 = F(x+dx3, y+dx3)[0]*dt
    dvy4 = F(x+dx3, y+dx3)[1]*dt
    dx4 = (vx+dvx3)*dt
    dy4 = (vy+dvy3)*dt
    

    # corrector
    vx += dvx1/6. + dvx2/3 + dvx2/3. + dvx4/6.
    vy += dvy1/6. + dvy2/3 + dvy2/3. + dvy4/6.
    x += dx1/6. + dx2/3. + dx3/3. + dx4/6.
    y += dy1/6. + dy2/3. + dy3/3. + dy4/6.
    
    t += dt
    
fgp.close()