Vlasiator ebf0dd394 on dev (v5.4.0 + 1054 commits)
Loading...
Searching...
No Matches
vectorpotentialdipole_fluxfunctions.py
Go to the documentation of this file.
1#!/usr/bin/env python import matplotlib.pyplot as plt
2
3# /*
4# * This file is part of Vlasiator.
5# * Copyright 2010-2016 Finnish Meteorological Institute
6# * Copyright 2017-2019 University of Helsinki
7# *
8# * For details of usage, see the COPYING file and read the "Rules of the Road"
9# * at http://www.physics.helsinki.fi/vlasiator/
10# *
11# * This program is free software; you can redistribute it and/or modify
12# * it under the terms of the GNU General Public License as published by
13# * the Free Software Foundation; either version 2 of the License, or
14# * (at your option) any later version.
15# *
16# * This program is distributed in the hope that it will be useful,
17# * but WITHOUT ANY WARRANTY; without even the implied warranty of
18# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19# * GNU General Public License for more details.
20# *
21# * You should have received a copy of the GNU General Public License along
22# * with this program; if not, write to the Free Software Foundation, Inc.,
23# * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24# */
25import numpy as np
26import math
27import sys,os
28import pytools as pt
29import matplotlib.pyplot as plt
30import matplotlib as mpl
31import scipy
32import fluxfunction as ff
33import fieldmodels
34
35''' Testing routine for different dipole formulations
36
37 Plots flux function contours of magnetic field in the meridional x-z-plane for four different models.
38 Also evaluates in-plane divergence, as the flux function doesn't work properly for the 3D dipole models.
39
40'''
41
42if len(sys.argv)!=1:
43 testset = int(sys.argv[1])
44else:
45 testset = 0
46
47plt.switch_backend('Agg')
48print(mpl.__version__)
49outfilename = "./vecpotdip_verify_fluxfunctions_"+str(testset)+".png"
50
51
52# Select flux function method
53ffc = ff.calculate # 3-way calculation
54#ffc = ff.calculate4 # 4-way calculation, returns mean of middle two values
55#ffc = ff.calculate4mean # 4-way calculation, returns mean of values
56
57#flux_levels = np.linspace(-5.e-5,5e-4,100)
58#flux_levels = np.linspace(-5.e-5,0,10000)
59
60#flux_levels = np.linspace(5.e-7,5e-4,1000)
61#flux_levels = np.linspace(5.e-12,5e-7,1000)
62#flux_levels = np.reshape([np.linspace(-1.e-5,-1e-14,200),np.linspace(1.e-14,1e-5,200)],400)
63#flux_levels = np.reshape([np.logspace(-14,-5,200)*-1,np.logspace(-14,-5,200)],400)
64
65#flux_levels = np.reshape([np.linspace(-1.e-6,-1e-14,300),np.linspace(1.e-14,1e-6,300)],600)
66flux_levels = np.reshape([np.linspace(-5.e-6,-1e-14,1000),np.linspace(1.e-14,5e-6,1000)],2000)
67
68RE=6371000.
69epsilon=1.e-15
70
71if testset==0:
72 tilt_angle_phi = 0.
73 tilt_angle_theta = 0.
74 BGB=[0,0,0]
75elif testset==1:
76 tilt_angle_phi = 10.
77 tilt_angle_theta = 0.
78 BGB=[0,0,0]
79elif testset==2:
80 tilt_angle_phi = 0.
81 tilt_angle_theta = 0.
82 BGB=[0,0,-5.e-9]
83elif testset==3: # Warning: this might not be a valid check, using flux functions with a out-of-plane tilted dipole?
84 tilt_angle_phi = 10.
85 tilt_angle_theta = 45.
86 BGB=[0,0,0]
87else: # Same as 0
88 print("Default")
89 tilt_angle_phi = 0.
90 tilt_angle_theta = 0.
91 BGB=[0,0,0]
92
93#fieldmodels.dipole.set_dipole(centerx, centery, centerz, tilt_phi, tilt_theta, mult=1.0, radius_f=None, radius_z=None):
94dip = fieldmodels.dipole(0,0,0,tilt_angle_phi,tilt_angle_theta)
95mdip = fieldmodels.dipole(80*RE,0,0,tilt_angle_phi,180.-tilt_angle_theta)
96
97
98fontsize=20
99# Create figure
100fig = plt.figure()
101fig.set_size_inches(20,20)
102
103gs = mpl.gridspec.GridSpec(2, 2, wspace=0.25, hspace=0.25)
104fig.add_subplot(gs[0, 0])
105fig.add_subplot(gs[0, 1])
106fig.add_subplot(gs[1, 0])
107fig.add_subplot(gs[1, 1])
108axes = fig.get_axes()
109
110fig.suptitle(r"Flux function contours of meridional plane magnetic field with dipole tilt $\Phi="+str(int(tilt_angle_phi))+"$, $\Theta="+str(int(tilt_angle_theta))+"$ with IMF=("+str(BGB[0])+","+str(BGB[1])+","+str(BGB[2])+")", fontsize=fontsize)
111
112
113nx = 200
114nz = 200
115xmin, xmax = (-59,41)
116zmin, zmax = (-50,50)
117
118x = np.linspace(xmin,xmax,num=nx)
119z = np.linspace(zmin,zmax,num=nz)
120BX = np.zeros([nx,1,nz])
121BZ = np.zeros([nx,1,nz])
122
123divB = np.zeros([nx,1,nz])
124
125[Xmesh,Zmesh] = scipy.meshgrid(x,z)
126
127dxdydz=[x[1]-x[0],0,z[1]-z[0]]
128
129
130ax = axes[0]
131print("0")
132for i in range(len(x)):
133 for j in range(len(z)):
134 BX[i,0,j] = dip.get_old(x[i]*RE,0,z[j]*RE,0,0,0) +BGB[0]
135 BZ[i,0,j] = dip.get_old(x[i]*RE,0,z[j]*RE,0,2,0) +BGB[2]
136 divB[i,0,j] = dip.get_old(x[i]*RE,0,z[j]*RE,1,0,0)
137 divB[i,0,j] += dip.get_old(x[i]*RE,0,z[j]*RE,1,2,2)
138print(np.sum(divB),np.amin(divB),np.amax(divB))
139ax.pcolormesh(Xmesh,Zmesh,divB[:,0,:])
140flux_function = ffc(BX,None,BZ,dxdydz)
141fluxcont = ax.contour(Xmesh,Zmesh,flux_function[:,0,:].T,flux_levels,colors='k',linestyles='solid',linewidths=0.5)
142ax.text(0.2,0.08,"Regular dipole",transform=ax.transAxes, bbox=dict(facecolor='white', alpha=0.7), fontsize=fontsize)
143
144ax = axes[1]
145print("1")
146for i in range(len(x)):
147 for j in range(len(z)):
148 BX[i,0,j] = dip.getX(x[i]*RE,0,z[j]*RE,0,0,0) +BGB[0]
149 BZ[i,0,j] = dip.getX(x[i]*RE,0,z[j]*RE,0,2,0) +BGB[2]
150 divB[i,0,j] = dip.getX(x[i]*RE,0,z[j]*RE,1,0,0)
151 divB[i,0,j] += dip.getX(x[i]*RE,0,z[j]*RE,1,2,2)
152print(np.sum(divB),np.amin(divB),np.amax(divB))
153ax.pcolormesh(Xmesh,Zmesh,divB[:,0,:])
154flux_function = ffc(BX,None,BZ,dxdydz)
155fluxcont = ax.contour(Xmesh,Zmesh,flux_function[:,0,:].T,flux_levels,colors='k',linestyles='solid',linewidths=0.5)
156ax.text(0.2,0.08,"Vector potential (X)",transform=ax.transAxes, bbox=dict(facecolor='white', alpha=0.7), fontsize=fontsize)
157
158# ax = axes[1]
159# print("1")
160# for i in range(len(x)):
161# for j in range(len(z)):
162# BX[i,0,j] = dip.getX(x[i]*RE,0,z[j]*RE,0,0,0)
163# BZ[i,0,j] = dip.getX(x[i]*RE,0,z[j]*RE,0,2,0)
164# BX[i,0,j] += mdip.getX(x[i]*RE,0,z[j]*RE,0,0,0) +BGB[0]
165# BZ[i,0,j] += mdip.getX(x[i]*RE,0,z[j]*RE,0,2,0) +BGB[2]
166# print(np.sum(divB),np.amin(divB),np.amax(divB))
167# ax.pcolormesh(Xmesh,Zmesh,divB[:,0,:])
168# flux_function = ffc(BX,None,BZ,dxdydz)
169# fluxcont = ax.contour(Xmesh,Zmesh,flux_function[:,0,:].T,flux_levels,colors='k',linestyles='solid',linewidths=0.5)
170# ax.text(0.2,0.08,"Vector potential + mirror (X)",transform=ax.transAxes, bbox=dict(facecolor='white', alpha=0.7), fontsize=fontsize)
171
172ax = axes[2]
173print("2")
174for i in range(len(x)):
175 for j in range(len(z)):
176 BX[i,0,j] = dip.get_old(x[i]*RE,0,z[j]*RE,0,0,0)
177 BZ[i,0,j] = dip.get_old(x[i]*RE,0,z[j]*RE,0,2,0)
178 BX[i,0,j] += mdip.get_old(x[i]*RE,0,z[j]*RE,0,0,0) +BGB[0]
179 BZ[i,0,j] += mdip.get_old(x[i]*RE,0,z[j]*RE,0,2,0) +BGB[2]
180
181 divB[i,0,j] = dip.get_old(x[i]*RE,0,z[j]*RE,1,0,0)
182 divB[i,0,j] += dip.get_old(x[i]*RE,0,z[j]*RE,1,2,2)
183 divB[i,0,j] += mdip.get_old(x[i]*RE,0,z[j]*RE,1,0,0)
184 divB[i,0,j] += mdip.get_old(x[i]*RE,0,z[j]*RE,1,2,2)
185print(np.sum(divB),np.amin(divB),np.amax(divB))
186ax.pcolormesh(Xmesh,Zmesh,divB[:,0,:])
187flux_function = ffc(BX,None,BZ,dxdydz)
188fluxcont = ax.contour(Xmesh,Zmesh,flux_function[:,0,:].T,flux_levels,colors='k',linestyles='solid',linewidths=0.5)
189ax.text(0.2,0.08,"Regular dipole + mirror",transform=ax.transAxes, bbox=dict(facecolor='white', alpha=0.7), fontsize=fontsize)
190
191if tilt_angle_phi<epsilon:
192 ax = axes[3]
193 print("3")
194 for i in range(len(x)):
195 for j in range(len(z)):
196 BX[i,0,j] = dip.get_ldp(x[i]*RE,0,z[j]*RE,0,0,0)
197 BZ[i,0,j] = dip.get_ldp(x[i]*RE,0,z[j]*RE,0,2,0)
198 BX[i,0,j] += mdip.get_ldp(x[i]*RE,0,z[j]*RE,0,0,0) +BGB[0]
199 BZ[i,0,j] += mdip.get_ldp(x[i]*RE,0,z[j]*RE,0,2,0) +BGB[2]
200
201 divB[i,0,j] = dip.get_ldp(x[i]*RE,0,z[j]*RE,1,0,0)
202 divB[i,0,j] += dip.get_ldp(x[i]*RE,0,z[j]*RE,1,2,2)
203 divB[i,0,j] += mdip.get_ldp(x[i]*RE,0,z[j]*RE,1,0,0)
204 divB[i,0,j] += mdip.get_ldp(x[i]*RE,0,z[j]*RE,1,2,2)
205 print(np.sum(divB),np.amin(divB),np.amax(divB))
206 ax.pcolormesh(Xmesh,Zmesh,divB[:,0,:])
207 flux_function = ffc(BX,None,BZ,dxdydz)
208 fluxcont = ax.contour(Xmesh,Zmesh,flux_function[:,0,:].T,flux_levels,colors='k',linestyles='solid',linewidths=0.5)
209 ax.text(0.2,0.08,"Line dipole + mirror",transform=ax.transAxes, bbox=dict(facecolor='white', alpha=0.7), fontsize=fontsize)
210
211fig.savefig(outfilename)
212plt.close()
vectorpotentialdipole_fluxfunctions.ffc
ffc
Definition vectorpotentialdipole_fluxfunctions.py:53