#!/usr/bin/env python
"""
Highly simplified model for the acceleration of the fast solar wind due to heating to a temperature T between x0 and max_r.

The derivation of the model is given as an exercise (P6_U4.pdf).

author: Bob Wimmer, IEAP/CAU Kiel, Germany, wimmer@physik.uni-kiel.de
date: June 21, 2012

"""



from numpy import *
from math import *
from scipy.integrate import odeint

kB = 1.38e-23
Msun = 1.99e30
Rsun = 6.96e8
G = 6.67e-11
amu=1.67e-27
mo = 16*amu
qe=1.6e-19
T=1.2e8
x0 = 2.*Rsun
max_r = 5.*Rsun


def func(y,dt):
    """This defines the acceleration in the corona
    """
    x,v = hsplit(y,2)
    xp = v
    if x < max_r:
        vp = 2.*kB*T/x/mo - G*Msun/x**2
    else:
        """Take into account adiabatic cooling of a sphere with radius gyroradius
        """
        vp = 2.*kB*T/x/mo*((x-Rsun)/(max_r-Rsun))**(-3./2.) - G*Msun/x**2
    #print x, vp
    z = hstack([xp, vp])
    return z


def main():
    outfile = open('speed_of_r.dat','w')
    vel = 0
    x   = x0
    i = 1.
    #adjust this to whatever you want. The i-bookkeeping is just for safety...
    while x < 2.*max_r and i < 100:
        t = [0.,i*10] #time step is 10s
        i +=1.
        y = hstack([x,vel])
        z = odeint(func,y,t)  #this is where the integration happens
        x,vel=hsplit(z[-1],2)
        string = str(x[0]/Rsun) + '\t' + str(vel[0]/1.e3) + '\n'
        outfile.write(string)
    outfile.close()
    
main()