I know I’m not very good, but according to your guidance and using the rest of the references and other guidance, I was able to change the code a little, maybe it’s a little difficult, but it gives me a good output.
I change the date in the date line.
According to the recorded number of sunspots, I ploted the maximum and minimum for 10 consecutive years and 36 consecutive months.
:Here I am writing the code although it is not very good
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from astropy import units as u
from astropy.coordinates import SkyCoord
from sunpy.map import Map
from sunpy.net import Fido, attrs as a
from sunkit_magex import pfss
import astropy.constants as const
class PFSSLoader:
def init(self, start_date, end_date, instrument=‘GONG’):
self.start_date = start_date
self.end_date = end_date
self.instrument = instrument
self.downloaded_files = None
self.gong_map = None
def load_data(self):
gong_fname = Fido.search(a.Time(self.start_date, self.end_date), a.Instrument(self.instrument))
self.downloaded_files = Fido.fetch(gong_fname)
self.gong_map = Map(self.downloaded_files[0])
def get_map(self):
return self.gong_map
class PFSSAnalyzer:
def init(self, gong_map, nrho=35, rss=2.5):
self.pfss_in = pfss.Input(gong_map, nrho, rss)
self.pfss_out = None
def compute_pfss(self):
self.pfss_out = pfss.pfss(self.pfss_in)
return self.pfss_out
def get_radial_field(self, ridx):
br = self.pfss_out.bc[0][:, :, ridx]
br_map = Map(br.T, self.pfss_out.source_surface_br.wcs)
r = np.exp(self.pfss_out.grid.rc[ridx])
return br_map, r
def trace_field_lines(self, r_factor=1.2, n_lat=8, n_lon=8):
tracer = pfss.tracing.FortranTracer()
r = r_factor * const.R_sun
lat = np.linspace(-np.pi / 2, np.pi / 2, n_lat, endpoint=False)
lon = np.linspace(0, 2 * np.pi, n_lon, endpoint=False)
lat, lon = np.meshgrid(lat, lon, indexing=‘ij’)
lat, lon = lat.ravel() * u.rad, lon.ravel() * u.rad
seeds = SkyCoord(lon, lat, r, frame=self.pfss_out.coordinate_frame)
return tracer.trace(seeds, self.pfss_out)
class PFSSPlotter:
@staticmethod
def set_axes_lims(ax):
ax.set_xlim(0, 360)
ax.set_ylim(0, 180)
def plot_input_field(self, pfss_input_map):
fig = plt.figure()
ax = plt.subplot(projection=pfss_input_map)
pfss_input_map.plot()
plt.colorbar()
ax.set_title(‘Input field’)
self.set_axes_lims(ax)
def plot_source_surface(self, pfss_output):
ss_br = pfss_output.source_surface_br
fig = plt.figure()
ax = plt.subplot(projection=ss_br)
ss_br.plot()
ax.plot_coord(pfss_output.source_surface_pils[0])
plt.colorbar()
ax.set_title(‘Source surface magnetic field’)
self.set_axes_lims(ax)
def plot_radial_field(self, br_map, r):
fig = plt.figure()
ax = plt.subplot(projection=br_map)
br_map.plot(cmap=‘RdBu’)
plt.colorbar()
ax.set_title(f’B_{{r}} at r={r:.2f}r_{{\\odot}}')
self.set_axes_lims(ax)
def plot_field_lines(self, field_lines):
fig = plt.figure()
ax = fig.add_subplot(111, projection=‘3d’)
for field_line in field_lines:
color = {0: ‘black’, -1: ‘tab:blue’, 1: ‘tab:red’}.get(field_line.polarity)
coords = field_line.coords
coords.representation_type = ‘cartesian’
ax.plot(coords.x / const.R_sun,
coords.y / const.R_sun,
coords.z / const.R_sun,
color=color, linewidth=1)
ax.set_title(‘PFSS solution’)
def show(self):
plt.show()
Example of usage:
if name == “main”:
Load data
start_date = input("Enter start date (YYYY/MM/DD HH:MM) | Example = 2021/12/08 00:00 : ")
end_date = input(“Enter end date (YYYY/MM/DD HH:MM) | Example = 2021/12/08 24:00 :”)
loader = PFSSLoader(start_date, end_date)
loader.load_data()
gong_map = loader.get_map()
Analyze data
analyzer = PFSSAnalyzer(gong_map)
pfss_out = analyzer.compute_pfss()
br_map, r = analyzer.get_radial_field(ridx=15)
field_lines = analyzer.trace_field_lines()
Plot data
plotter = PFSSPlotter()
plotter.plot_input_field(analyzer.pfss_in.map)
plotter.plot_source_surface(pfss_out)
plotter.plot_radial_field(br_map, r)
plotter.plot_field_lines(field_lines)
plotter.show()
Nabil Freij via OpenAstronomy Community <notifications@openastronomy.discoursemail.com> در تاریخ چهارشنبه ۲۸ اوت ۲۰۲۴ ساعت ۱۶:۳۵ نوشت: