diff --git a/.github/workflows/test_main.yml b/.github/workflows/test_main.yml new file mode 100644 index 0000000..d2803a5 --- /dev/null +++ b/.github/workflows/test_main.yml @@ -0,0 +1,46 @@ +name: Python package + +on: [push] + +jobs: + build: + + runs-on: ubuntu-latest + strategy: + max-parallel: 4 + matrix: + python-version: [3.6, 3.7, 3.8, 3.9] + + steps: + - uses: actions/checkout@v1 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v1 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -e . + pip install -r requirements_dev.txt + - name: Lint with flake8 + run: | + pip install flake8 + # stop the build if there are Python syntax errors or undefined names + flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics + # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide + flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics + - name: Test with pytest + run: | + pip install pytest + pip install pytest-cov + pytest --cov=./ --cov-report=xml + - name: Upload coverage to Codecov + if: matrix.python-version == '3.8' + uses: codecov/codecov-action@v1 + with: + token: ${{ secrets.CODECOV_TOKEN }} + file: ./coverage.xml + flags: unittests + name: codecov-umbrella + yml: ./codecov.yml + fail_ci_if_error: true diff --git a/.github/workflows/test_pull_requests.yml b/.github/workflows/test_pull_requests.yml new file mode 100644 index 0000000..b14cc41 --- /dev/null +++ b/.github/workflows/test_pull_requests.yml @@ -0,0 +1,46 @@ +name: Python package + +on: [pull_request] + +jobs: + build: + + runs-on: ubuntu-latest + strategy: + max-parallel: 4 + matrix: + python-version: [3.6, 3.7, 3.8, 3.9] + + steps: + - uses: actions/checkout@v1 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v1 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -e . + pip install -r requirements_dev.txt + - name: Lint with flake8 + run: | + pip install flake8 + # stop the build if there are Python syntax errors or undefined names + flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics + # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide + flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics + - name: Test with pytest + run: | + pip install pytest + pip install pytest-cov + pytest --cov=./ --cov-report=xml + - name: Upload coverage to Codecov + if: matrix.python-version == '3.8' + uses: codecov/codecov-action@v1 + with: + token: ${{ secrets.CODECOV_TOKEN }} + file: ./coverage.xml + flags: unittests + name: codecov-umbrella + yml: ./codecov.yml + fail_ci_if_error: true diff --git a/README.rst b/README.rst index f8d3c8c..0fa27a4 100644 --- a/README.rst +++ b/README.rst @@ -13,6 +13,9 @@ space :target: https://space.readthedocs.io/en/latest/?badge=latest :alt: Documentation Status +.. image:: https://codecov.io/gh/LaboratoryOfPlasmaPhysics/space/branch/main/graph/badge.svg?branch=main + :target: https://codecov.io/gh/LaboratoryOfPlasmaPhysics/space/branch/main + :alt: Coverage Status diff --git a/requirements_dev.txt b/requirements_dev.txt index c54c8ce..72c2ca6 100644 --- a/requirements_dev.txt +++ b/requirements_dev.txt @@ -8,3 +8,7 @@ coverage==4.5.4 Sphinx==1.8.5 twine==1.14.0 +ddt +pandas +numpy +matplotlib diff --git a/setup.py b/setup.py index a59c67d..a2777fb 100644 --- a/setup.py +++ b/setup.py @@ -10,7 +10,7 @@ with open('HISTORY.rst') as history_file: history = history_file.read() -requirements = ['pandas', 'numpy', 'matplotlib'] +requirements = ['pandas', 'numpy', 'matplotlib', 'scipy'] setup_requirements = [] diff --git a/space/coordinates/coordinates.py b/space/coordinates/coordinates.py index ac1bed3..124f905 100644 --- a/space/coordinates/coordinates.py +++ b/space/coordinates/coordinates.py @@ -1,91 +1,100 @@ import numpy as np import pandas as pd -from datetime import timedelta, datetime + def spherical_to_cartesian(R, theta, phi): x = R * np.cos(theta) - y = R * np.sin(theta)*np.cos(phi) - z = R * np.sin(theta)*np.sin(phi) + y = R * np.sin(theta) * np.cos(phi) + z = R * np.sin(theta) * np.sin(phi) return x, y, z -def cartesian_to_spherical(X,Y,Z): - r = np.sqrt(X**2+Y**2+Z**2) - theta = np.arccos(X/r) - phi = np.arctan2(Z,Y) - return r,theta,phi - -def BaseChoice(base,R,theta,phi): - if base=='cartesian' : - return spherical_to_cartesian(R,theta,phi) - elif base=='spherical' : - return R,theta,phi - else : + +def cartesian_to_spherical(X, Y, Z): + r = np.sqrt(X ** 2 + Y ** 2 + Z ** 2) + theta = np.arccos(X / r) + phi = np.arctan2(Z, Y) + return r, theta, phi + + +def base_choice(base, R, theta, phi): + if base == 'cartesian': + return spherical_to_cartesian(R, theta, phi) + elif base == 'spherical': + return R, theta, phi + else: print('Error : base parameter must be set to "cartesian" or "spherical" ') -def Check_base(pos): - if ((hasattr(pos, 'R' )) | (hasattr(pos, 'r' ))) & (hasattr(pos, 'X' )) | (hasattr(pos, 'x' )): +def check_base(pos): + if ((hasattr(pos, 'R')) | (hasattr(pos, 'r'))) & (hasattr(pos, 'X')) | (hasattr(pos, 'x')): base = 'spherical&cartesian' - elif (hasattr(pos, 'R' )) | (hasattr(pos, 'r' )) : + elif (hasattr(pos, 'R')) | (hasattr(pos, 'r')): base = 'spherical' - elif (hasattr(pos, 'X' )) | (hasattr(pos, 'x' )) : + elif (hasattr(pos, 'X')) | (hasattr(pos, 'x')): base = 'cartesian' - else : + else: print('must check the name of the variables : cartesian = (X,Y,Z) and spherical = (R,theta,phi)') return base -def Add_cst_radiuus(pos,cst): - base = Check_base(pos) - if base=='cartesian': - r,theta,phi = CartesianToSpherical(pos.X,pos.Y,pos.Z) - else : - r,theta,phi = pos.R, pos.theta, pos.phi - r=r+cst - x,y,z = spherical_to_cartesian(r,theta,phi) - return pd.DataFrame({'X' : x, 'Y' : y, 'Z' :z}) -def SWI_base(omni_data): - norm_V =np.sqrt(omni_data.Vx**2+(omni_data.Vy+29.8)**2+omni_data.Vz**2) - X = np.array([-np.array([vx,vy+29.8,vz])/V for vx,vy,vz,V in zip( omni_data.Vx.values, omni_data.Vy.values, omni_data.Vz.values,norm_V)]) - B = np.array([np.array([bx,by,bz]) for bx,by,bz in zip( np.sign(omni_data.COA.values)*omni_data.Bx.values, np.sign(omni_data.COA.values)*omni_data.By.values, np.sign(omni_data.COA.values)*omni_data.Bz.values)]) +def add_cst_radiuus(pos, cst): + base = check_base(pos) + if base == 'cartesian': + r, theta, phi = cartesian_to_spherical(pos.X, pos.Y, pos.Z) + else: + r, theta, phi = pos.R, pos.theta, pos.phi + r = r + cst + x, y, z = spherical_to_cartesian(r, theta, phi) + return pd.DataFrame({'X': x, 'Y': y, 'Z': z}) - Z = np.cross(X,B) +def swi_base(omni_data): + norm_V = np.sqrt(omni_data.Vx ** 2 + (omni_data.Vy + 29.8) ** 2 + omni_data.Vz ** 2) + X = np.array([-np.array([vx, vy + 29.8, vz]) / V for vx, vy, vz, V in + zip(omni_data.Vx.values, omni_data.Vy.values, omni_data.Vz.values, norm_V)]) + B = np.array([np.array([bx, by, bz]) for bx, by, bz in zip(np.sign(omni_data.COA.values) * omni_data.Bx.values, + np.sign(omni_data.COA.values) * omni_data.By.values, + np.sign(omni_data.COA.values) * omni_data.Bz.values)]) - norm_cross = np.array(np.sqrt(Z[:,0]**2+Z[:,1]**2+Z[:,2]**2)) - Z = np.array([z/n for z,n in zip(Z,norm_cross)]) - Y = np.cross(Z,X) - return X,Y,Z + Z = np.cross(X, B) + + norm_cross = np.array(np.sqrt(Z[:, 0] ** 2 + Z[:, 1] ** 2 + Z[:, 2] ** 2)) + Z = np.array([z / n for z, n in zip(Z, norm_cross)]) + Y = np.cross(Z, X) + return X, Y, Z def to_swi(omni_data, msh_data, pos_msh): - X_swi, Y_swi, Z_swi = SWI_base(omni_data) + X_swi, Y_swi, Z_swi = swi_base(omni_data) data = msh_data.copy() o_data = omni_data.copy() pos = pos_msh.copy() - - o_data['Vx'] = X_swi[:,0]*omni_data['Vx']+X_swi[:,1]*omni_data['Vy']+X_swi[:,2]*omni_data['Vz'] - o_data['Vy'] = Y_swi[:,0]*omni_data['Vx']+Y_swi[:,1]*omni_data['Vy']+Y_swi[:,2]*omni_data['Vz'] - o_data['Vz'] = Z_swi[:,0]*omni_data['Vx']+Z_swi[:,1]*omni_data['Vy']+Z_swi[:,2]*omni_data['Vz'] - - o_data['Bx'] = np.sign(omni_data['COA'])*(X_swi[:,0]*omni_data['Bx']+X_swi[:,1]*omni_data['By']+X_swi[:,2]*omni_data['Bz']) - o_data['By'] = np.sign(omni_data['COA'])*(Y_swi[:,0]*omni_data['Bx']+Y_swi[:,1]*omni_data['By']+Y_swi[:,2]*omni_data['Bz']) - o_data['Bz'] = np.sign(omni_data['COA'])*(Z_swi[:,0]*omni_data['Bx']+Z_swi[:,1]*omni_data['By']+Z_swi[:,2]*omni_data['Bz']) - - data['Vx'] = X_swi[:,0]*msh_data['Vx']+X_swi[:,1]*msh_data['Vy']+X_swi[:,2]*msh_data['Vz'] - data['Vy'] = Y_swi[:,0]*msh_data['Vx']+Y_swi[:,1]*msh_data['Vy']+Y_swi[:,2]*msh_data['Vz'] - data['Vz'] = Z_swi[:,0]*msh_data['Vx']+Z_swi[:,1]*msh_data['Vy']+Z_swi[:,2]*msh_data['Vz'] - - data['Bx'] = np.sign(omni_data['COA'])*(X_swi[:,0]*msh_data['Bx']+X_swi[:,1]*msh_data['By']+X_swi[:,2]*msh_data['Bz']) - data['By'] = np.sign(omni_data['COA'])*(Y_swi[:,0]*msh_data['Bx']+Y_swi[:,1]*msh_data['By']+Y_swi[:,2]*msh_data['Bz']) - data['Bz'] = np.sign(omni_data['COA'])*(Z_swi[:,0]*msh_data['Bx']+Z_swi[:,1]*msh_data['By']+Z_swi[:,2]*msh_data['Bz']) - - - pos['X'] = X_swi[:,0]*pos_msh['X']+X_swi[:,1]*pos_msh['Y']+X_swi[:,2]*pos_msh['Z'] - pos['Y'] = Y_swi[:,0]*pos_msh['X']+Y_swi[:,1]*pos_msh['Y']+Y_swi[:,2]*pos_msh['Z'] - pos['Z'] = Z_swi[:,0]*pos_msh['X']+Z_swi[:,1]*pos_msh['Y']+Z_swi[:,2]*pos_msh['Z'] - - - return data,pos,o_data + o_data['Vx'] = X_swi[:, 0] * omni_data['Vx'] + X_swi[:, 1] * omni_data['Vy'] + X_swi[:, 2] * omni_data['Vz'] + o_data['Vy'] = Y_swi[:, 0] * omni_data['Vx'] + Y_swi[:, 1] * omni_data['Vy'] + Y_swi[:, 2] * omni_data['Vz'] + o_data['Vz'] = Z_swi[:, 0] * omni_data['Vx'] + Z_swi[:, 1] * omni_data['Vy'] + Z_swi[:, 2] * omni_data['Vz'] + + o_data['Bx'] = np.sign(omni_data['COA']) * ( + X_swi[:, 0] * omni_data['Bx'] + X_swi[:, 1] * omni_data['By'] + X_swi[:, 2] * omni_data['Bz']) + o_data['By'] = np.sign(omni_data['COA']) * ( + Y_swi[:, 0] * omni_data['Bx'] + Y_swi[:, 1] * omni_data['By'] + Y_swi[:, 2] * omni_data['Bz']) + o_data['Bz'] = np.sign(omni_data['COA']) * ( + Z_swi[:, 0] * omni_data['Bx'] + Z_swi[:, 1] * omni_data['By'] + Z_swi[:, 2] * omni_data['Bz']) + + data['Vx'] = X_swi[:, 0] * msh_data['Vx'] + X_swi[:, 1] * msh_data['Vy'] + X_swi[:, 2] * msh_data['Vz'] + data['Vy'] = Y_swi[:, 0] * msh_data['Vx'] + Y_swi[:, 1] * msh_data['Vy'] + Y_swi[:, 2] * msh_data['Vz'] + data['Vz'] = Z_swi[:, 0] * msh_data['Vx'] + Z_swi[:, 1] * msh_data['Vy'] + Z_swi[:, 2] * msh_data['Vz'] + + data['Bx'] = np.sign(omni_data['COA']) * ( + X_swi[:, 0] * msh_data['Bx'] + X_swi[:, 1] * msh_data['By'] + X_swi[:, 2] * msh_data['Bz']) + data['By'] = np.sign(omni_data['COA']) * ( + Y_swi[:, 0] * msh_data['Bx'] + Y_swi[:, 1] * msh_data['By'] + Y_swi[:, 2] * msh_data['Bz']) + data['Bz'] = np.sign(omni_data['COA']) * ( + Z_swi[:, 0] * msh_data['Bx'] + Z_swi[:, 1] * msh_data['By'] + Z_swi[:, 2] * msh_data['Bz']) + + pos['X'] = X_swi[:, 0] * pos_msh['X'] + X_swi[:, 1] * pos_msh['Y'] + X_swi[:, 2] * pos_msh['Z'] + pos['Y'] = Y_swi[:, 0] * pos_msh['X'] + Y_swi[:, 1] * pos_msh['Y'] + Y_swi[:, 2] * pos_msh['Z'] + pos['Z'] = Z_swi[:, 0] * pos_msh['X'] + Z_swi[:, 1] * pos_msh['Y'] + Z_swi[:, 2] * pos_msh['Z'] + + return data, pos, o_data diff --git a/space/models/planetary.py b/space/models/planetary.py index 95027a6..8943366 100644 --- a/space/models/planetary.py +++ b/space/models/planetary.py @@ -3,32 +3,30 @@ from scipy import constants as cst import sys + sys.path.append('.') -from .. import utils from ..coordinates import coordinates as coords from ..smath import resolve_poly2 - def _checking_angles(theta, phi): - if isinstance(theta, np.ndarray) and isinstance(theta, np.ndarray) and len(theta.shape) > 1 and len(phi.shape) > 1: return np.meshgrid(theta, phi) return theta, phi def _formisano1979(theta, phi, **kwargs): - a11,a22,a33,a12,a13,a23,a14,a24,a34,a44 = kwargs["coefs"] + a11, a22, a33, a12, a13, a23, a14, a24, a34, a44 = kwargs["coefs"] ct = np.cos(theta) st = np.sin(theta) cp = np.cos(phi) sp = np.sin(phi) - A = a11 * ct**2 + st**2 * (a22* cp**2 + a33 * sp**2) \ - + ct * st * (a12 * cp + a13 * sp) + a23 * st**2 * cp * sp - B = a14*ct + st*(a24*cp + a34*sp) + A = a11 * ct ** 2 + st ** 2 * (a22 * cp ** 2 + a33 * sp ** 2) \ + + ct * st * (a12 * cp + a13 * sp) + a23 * st ** 2 * cp * sp + B = a14 * ct + st * (a24 * cp + a34 * sp) C = a44 - D = B**2 - 4*A*C - return (-B + np.sqrt(D))/(2*A) + D = B ** 2 - 4 * A * C + return (-B + np.sqrt(D)) / (2 * A) def formisano1979(theta, phi, **kwargs): @@ -52,15 +50,15 @@ def formisano1979(theta, phi, **kwargs): ''' if kwargs["boundary"] == "magnetopause": - coefs = [0.65,1,1.16,0.03,-0.28,-0.11,21.41,0.46,-0.36,-221] + coefs = [0.65, 1, 1.16, 0.03, -0.28, -0.11, 21.41, 0.46, -0.36, -221] elif kwargs["boundary"] == "bow_shock": coefs = [0.52, 1, 1.05, 0.13, -0.16, -0.08, 47.53, -0.42, 0.67, -613] else: raise ValueError("boundary: {} not allowed".format(kwargs["boundary"])) theta, phi = _checking_angles(theta, phi) - r = _formisano1979(theta, phi, coefs = coefs) - base = kwargs.get("base", "cartesian") + r = _formisano1979(theta, phi, coefs=coefs) + base = kwargs.get("base", "cartesian") if base == "cartesian": return coords.spherical_to_cartesian(r, theta, phi) elif base == "spherical": @@ -68,18 +66,15 @@ def formisano1979(theta, phi, **kwargs): raise ValueError("unknown base '{}'".format(kwargs["base"])) - def mp_formisano1979(theta, phi, **kwargs): return formisano1979(theta, phi, boundary="magnetopause", **kwargs) + def bs_formisano1979(theta, phi, **kwargs): return formisano1979(theta, phi, boundary="bow_shock", **kwargs) - - def Fairfield1971(x, args): - ''' Fairfield 1971 : Magnetopause and Bow shock models. Give positions of the boudaries in plans (XY) with Z=0 and (XZ) with Y=0. function's arguments : @@ -93,21 +88,20 @@ def Fairfield1971(x, args): return : DataFrame (Pandas) with the position (X,Y,Z) in Re of the wanted boudary to plot (XY) and (XZ) plans. ''' - - A,B,C,D,E = args[0],args[1],args[2],args[3],args[4] + A, B, C, D, E = args[0], args[1], args[2], args[3], args[4] a = 1 - b = A*x + C - c = B*x**2 + D*x + E + b = A * x + C + c = B * x ** 2 + D * x + E - delta = b**2-4*a*c + delta = b ** 2 - 4 * a * c - ym = (-b - np.sqrt(delta))/(2*a) - yp = (-b + np.sqrt(delta))/(2*a) + ym = (-b - np.sqrt(delta)) / (2 * a) + yp = (-b + np.sqrt(delta)) / (2 * a) - pos=pd.DataFrame({'X' : np.concatenate([x, x[::-1]]), - 'Y' : np.concatenate([yp, ym[::-1]]), - 'Z' : np.concatenate([yp, ym[::-1]]),}) + pos = pd.DataFrame({'X': np.concatenate([x, x[::-1]]), + 'Y': np.concatenate([yp, ym[::-1]]), + 'Z': np.concatenate([yp, ym[::-1]]), }) return pos.dropna() @@ -156,19 +150,16 @@ def make_Rav(theta, phi): D = 0.937 * (0.846 + 0.042 * B) R0 = make_Rav(0, 0) - - Rav = make_Rav(theta, phi) K = ((gamma - 1) * Ma ** 2 + 2) / ((gamma + 1) * (Ma ** 2 - 1)) r = (Rav / R0) * (C / (Np * V ** 2) ** (1 / 6)) * (1 + D * K) - base = kwargs.get('base', 'cartesian') if base == "cartesian": x = r * np.cos(theta) - y = r * np.sin(theta)*np.cos(phi) - z = r * np.sin(theta)*np.sin(phi) + y = r * np.sin(theta) * np.cos(phi) + z = r * np.sin(theta) * np.sin(phi) return x, y, z elif base == "spherical": return r, theta, phi @@ -195,10 +186,9 @@ def mp_shue1998(theta, phi, **kwargs): Pd = kwargs.get("pdyn", 2) Bz = kwargs.get("Bz", 1) - r0 = (10.22 + 1.29 * np.tanh(0.184 * (Bz + 8.14))) * Pd**(-1./6.6) - a = (0.58-0.007*Bz)*(1+0.024*np.log(Pd)) - r = r0*(2./(1+np.cos(theta)))**a - + r0 = (10.22 + 1.29 * np.tanh(0.184 * (Bz + 8.14))) * Pd ** (-1. / 6.6) + a = (0.58 - 0.007 * Bz) * (1 + 0.024 * np.log(Pd)) + r = r0 * (2. / (1 + np.cos(theta))) ** a base = kwargs.get("base", "cartesian") if base == "cartesian": @@ -211,8 +201,7 @@ def mp_shue1998(theta, phi, **kwargs): raise ValueError("unknown base '{}'".format(kwargs["base"])) - -def MP_Lin2010(phi_in ,th_in, Pd, Pm, Bz, tilt=0.): +def MP_Lin2010(phi_in, th_in, Pd, Pm, Bz, tilt=0.): ''' The Lin 2010 Magnetopause model. Returns the MP distance for a given azimuth (phi), zenith (th), solar wind dynamic and magnetic pressures (nPa) and Bz (in nT). @@ -247,48 +236,48 @@ def MP_Lin2010(phi_in ,th_in, Pd, Pm, Bz, tilt=0.): arr = type(np.array([])) - if(type(th_in) == arr): + if (type(th_in) == arr): th = th_in.copy() else: th = th_in - if(type(phi_in) == arr): + if (type(phi_in) == arr): phi = phi_in.copy() else: phi = phi_in el = th_in < 0. - if(type(el) == arr): - if(el.any()): + if (type(el) == arr): + if (el.any()): th[el] = -th[el] - if(type(phi) == type(arr)): - phi[el] = phi[el]+np.pi + if (type(phi) == type(arr)): + phi[el] = phi[el] + np.pi else: - phi = phi*np.ones(th.shape)+np.pi*el + phi = phi * np.ones(th.shape) + np.pi * el else: - if(el): + if (el): th = -th - phi = phi+np.pi + phi = phi + np.pi - P = Pd+Pm + P = Pd + Pm def exp2(i): - return a[i]*(np.exp(a[i+1]*Bz)-1)/(np.exp(a[i+2]*Bz)+1) + return a[i] * (np.exp(a[i + 1] * Bz) - 1) / (np.exp(a[i + 2] * Bz) + 1) def quad(i, s): - return a[i]+s[0]*a[i+1]*tilt+s[1]*a[i+2]*tilt**2 + return a[i] + s[0] * a[i + 1] * tilt + s[1] * a[i + 2] * tilt ** 2 - r0 = a[0]*P**a[1]*(1+exp2(2)) + r0 = a[0] * P ** a[1] * (1 + exp2(2)) beta = [a[6] + exp2(7), a[10], quad(11, [1, 0]), a[13]] - f = np.cos(0.5*th)+a[5]*np.sin(2*th)*(1-np.exp(-th)) - s = beta[0]+beta[1]*np.sin(phi)+beta[2]*np.cos(phi)+beta[3]*np.cos(phi)**2 - f = f**(s) + f = np.cos(0.5 * th) + a[5] * np.sin(2 * th) * (1 - np.exp(-th)) + s = beta[0] + beta[1] * np.sin(phi) + beta[2] * np.cos(phi) + beta[3] * np.cos(phi) ** 2 + f = f ** (s) c = {} d = {} @@ -296,19 +285,20 @@ def quad(i, s): PHI = {} e = {} for i, s in zip(['n', 's'], [1, -1]): - c[i] = a[14]*P**a[15] + c[i] = a[14] * P ** a[15] d[i] = quad(16, [s, 1]) TH[i] = quad(19, [s, 0]) - PHI[i] = np.cos(th)*np.cos(TH[i]) - PHI[i] = PHI[i] + np.sin(th)*np.sin(TH[i])*np.cos(phi-(1-s)*0.5*np.pi) + PHI[i] = np.cos(th) * np.cos(TH[i]) + PHI[i] = PHI[i] + np.sin(th) * np.sin(TH[i]) * np.cos(phi - (1 - s) * 0.5 * np.pi) PHI[i] = np.arccos(PHI[i]) e[i] = a[21] - r = f*r0 + r = f * r0 + + Q = c['n'] * np.exp(d['n'] * PHI['n'] ** e['n']) + Q = Q + c['s'] * np.exp(d['s'] * PHI['s'] ** e['s']) - Q = c['n']*np.exp(d['n']*PHI['n']**e['n']) - Q = Q + c['s']*np.exp(d['s']*PHI['s']**e['s']) + return r + Q - return r+Q def mp_liu2015(theta, phi, **kwargs): if isinstance(theta, np.ndarray) | isinstance(theta, pd.Series): @@ -384,11 +374,9 @@ def mp_liu2015(theta, phi, **kwargs): raise ValueError("unknown base '{}'".format(kwargs["base"])) - - _models = {"mp_shue1998": mp_shue1998, "mp_formisano1979": mp_formisano1979, - "mp_liu2015" : mp_liu2015, + "mp_liu2015": mp_liu2015, "bs_formisano1979": bs_formisano1979, "bs_jerab2005": bs_Jerab2005} @@ -417,8 +405,7 @@ def _parabolic_approx(theta, phi, x, xf, **kwargs): b = 4 * K * np.cos(theta) c = -4 * K * x r = resolve_poly2(a, b, c, 0) - return coords.BaseChoice(kwargs.get("base", "cartesian"), r, theta, phi) - + return coords.base_choice(kwargs.get("base", "cartesian"), r, theta, phi) def check_parabconfoc(func): @@ -428,6 +415,7 @@ def wrapper(self, theta, phi, **kwargs): if kwargs["parabolic"] is False and kwargs["confocal"] is True: raise ValueError("cannot be confocal if not parabolic") return func(self, theta, phi, **kwargs) + return wrapper @@ -453,7 +441,6 @@ def magnetopause(self, theta, phi, **kwargs): else: return self._magnetopause(theta, phi, **kwargs) - @check_parabconfoc def bow_shock(self, theta, phi, **kwargs): if kwargs["parabolic"]: @@ -466,10 +453,9 @@ def boundaries(self, theta, phi, **kwargs): if kwargs["parabolic"]: return self._parabolize(theta, phi, **kwargs) else: - return self._magnetopause(theta, phi, **kwargs),\ + return self._magnetopause(theta, phi, **kwargs), \ self._bow_shock(theta, phi, **kwargs) - def _parabolize(self, theta, phi, **kwargs): xmp, y, xfmp = _interest_points(self._magnetopause, **kwargs) xbs, y, xfbs = _interest_points(self._bow_shock, **kwargs) diff --git a/space/plot/planet_env.py b/space/plot/planet_env.py index 5810ffb..026494d 100644 --- a/space/plot/planet_env.py +++ b/space/plot/planet_env.py @@ -1,12 +1,18 @@ - import matplotlib.pyplot as plt +import numpy as np + def layout_EarthEnv_3planes(**kwargs): - #figsize=(15,4.5),xlim=(-30,30),ylim=(-30,30),zlim =(-30,30),alpha=0.5): + # figsize=(15,4.5), + # to be fixed! + xlim = (-30, 30) + ylim = (-30, 30) + zlim = (-30, 30) + alpha = 0.5 - figsize = kwargs.get("kwargs", (15, 4.5)) + figsize = kwargs.get("kwargs", (15, 4.5)) # kwargs?! - fig, (ax0, ax1, ax2) = plt.subplots(ncols=3,figsize=figsize,constrained_layout=True) + fig, (ax0, ax1, ax2) = plt.subplots(ncols=3, figsize=figsize, constrained_layout=True) ax0.set_xlabel('X (Re)') ax0.set_ylabel('Y (Re)') @@ -27,33 +33,31 @@ def layout_EarthEnv_3planes(**kwargs): ax2.set_xlim(ylim) ax2.set_ylim(zlim) - return ax0,ax1,ax2 + return ax0, ax1, ax2 + # plot_boundaries(MP, BS, slice_x=22, slice_y=24, slice_z=0) def make_YZ_plan(pos): - a = np.linspace(0,2*np.pi,100) - r=abs(pos[(pos.X**2).argmin(): (pos.X**2).argmin()+1].Y.values) - return(r*np.cos(a),r*np.sin(a)) - - - + a = np.linspace(0, 2 * np.pi, 100) + r = abs(pos[(pos.X ** 2).argmin(): (pos.X ** 2).argmin() + 1].Y.values) + return r * np.cos(a), r * np.sin(a) -def plot_boudaries(MP, BS, **kwargs): - style = kwargs.get("style", ['--k' , '--k']) +def plot_boundaries(MP, BS, **kwargs): + style = kwargs.get("style", ['--k', '--k']) alpha = kwargs.get("alpha", 0.6) axes = kwargs.get("axes", layout_EarthEnv_3planes(**kwargs)) if "slice_x" in kwargs: axes[0].plot() - axes[0].plot(MP.X,MP.Y,style[0], alpha=alpha) - axes[0].plot(BS.X,BS.Y,style[1], alpha=alpha) + axes[0].plot(MP.X, MP.Y, style[0], alpha=alpha) + axes[0].plot(BS.X, BS.Y, style[1], alpha=alpha) - axes[1].plot(MP.X,MP.Z,style[0], alpha=alpha) - axes[1].plot(BS.X,BS.Z,style[1], alpha=alpha) + axes[1].plot(MP.X, MP.Z, style[0], alpha=alpha) + axes[1].plot(BS.X, BS.Z, style[1], alpha=alpha) - axes[2].plot(make_YZ_plan(MP)[0], make_YZ_plan(MP)[1],style[0],alpha=alpha) - axes[2].plot(make_YZ_plan(BS)[0], make_YZ_plan(BS)[1],style[1],alpha=alpha) + axes[2].plot(make_YZ_plan(MP)[0], make_YZ_plan(MP)[1], style[0], alpha=alpha) + axes[2].plot(make_YZ_plan(BS)[0], make_YZ_plan(BS)[1], style[1], alpha=alpha)