Commit e4f879c8 authored by Qusai Al Shidi's avatar Qusai Al Shidi 💬
Browse files

Transformation function made. Write IMF untested. Almost done.

parent 0d243de1
......@@ -6,10 +6,11 @@ __author__ = 'Qusai Al Shidi'
__email__ = 'qusai@umich.edu'
import numpy as np
from scipy.optimize import newton, check_grad
def get_v_gsm_from_omni(dict_omni):
"""Adds gsm values to the omni dictionary
"""Adds gsm values to the omni dictionary
This is helpful for the omni dict created with #swmfpy.web.get_omni_data()
......@@ -40,13 +41,76 @@ def get_v_gsm_from_omni(dict_omni):
# Author: Qusai Al Shidi
# Email: qusai@umich.edu
vy_gse, vz_gse = dict_omni['vy_gse'], dict_omni['vz_gse']
by_gse, bz_gse = dict_omni['by_gse'], dict_omni['bz_gse']
by_gsm, bz_gsm = dict_omni['by'], dict_omni['bz']
velocity = np.sqrt(vx_gse**2 + vz_gse**2)
theta = np.arctan(bz_gsm, by_gsm)
dict_omni['vy'] = velocity*np.cos(theta)
dict_omni['vz'] = velocity*np.sin(theta)
b_gsm = [dict_omni['bx'], dict_omni['by'], dict_omni['bz']]
b_gse = [dict_omni['bx'], dict_omni['by_gse'], dict_omni['bz_gse']]
v_gse = [dict_omni['vx_gse'], dict_omni['vy_gse'], dict_omni['vz_gse']]
def check_prime(func, fprime, *args):
'Make sure my derivatives are correct'
assert check_grad(func, fprime, [0.1], *args) < 1.e-3, 'bad grad'
def rot_y(theta, vec):
'rotate y in yz plane'
return vec[1]*np.cos(theta)-vec[2]*np.sin(theta)
def rot_y_prime(theta, vec):
'rotate y in yz plane'
return -vec[1]*np.sin(theta)-vec[2]*np.cos(theta)
check_prime(rot_y, rot_y_prime, (1, 1, 1))
def rot_y_prime2(theta, vec):
'rotate y in yz plane'
return -vec[1]*np.cos(theta)+vec[2]*np.sin(theta)
check_prime(rot_y_prime, rot_y_prime2, (1, 1, 1))
def rot_z(theta, vec):
'rotate z in yz plane'
return vec[1]*np.sin(theta)+vec[2]*np.cos(theta)
def rot_z_prime(theta, vec):
'rotate z in yz plane'
return vec[1]*np.cos(theta)-vec[2]*np.sin(theta)
check_prime(rot_z, rot_z_prime, (1, 1, 1))
def rot_z_prime2(theta, vec):
'rotate z in yz plane'
return -vec[1]*np.sin(theta)-vec[2]*np.cos(theta)
check_prime(rot_z_prime, rot_z_prime2, (1, 1, 1))
def opt_func(theta, vec_p, vec):
'The two functions to optimize'
opt_y = vec_p[1] - rot_y(theta, vec)
opt_z = vec_p[2] - rot_z(theta, vec)
return opt_y + opt_z
def opt_func_prime(theta, _, vec):
'The two functions to optimize'
opt_y = -rot_y_prime(theta, vec)
opt_z = -rot_z_prime(theta, vec)
return opt_y + opt_z
check_prime(opt_func, opt_func_prime, (1, 1, 1), (1, 1, 1))
def opt_func_prime2(theta, _, vec):
'The two functions to optimize'
opt_y = -rot_y_prime2(theta, vec)
opt_z = -rot_z_prime2(theta, vec)
return opt_y + opt_z
check_prime(opt_func_prime, opt_func_prime2, (1, 1, 1), (1, 1, 1))
angle = newton(opt_func, np.zeros_like(b_gsm[0]), args=(b_gsm, b_gse),
fprime=opt_func_prime, fprime2=opt_func_prime2, tol=1.e-2)
dict_omni.update({
'vx': dict_omni['vx_gse'],
'vy': rot_y(angle, v_gse),
'vz': rot_z(angle, v_gse),
})
return dict_omni
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment