I still haven’t been able to plot such a model and my codes completely failed in Python
Hello @Bahar,
Sorry to hear that it keeps failing. I can’t help you as I have no information on what you are using in Python and what is failing. If you could provide that information, we can go from there.
In the meantime, there are a set of examples here: Examples — sunkit-magex 1.0.0 documentation
It might be possible that there something there that would be useful for you.
Thanks for the tip, I’ve used the solved examples before and they worked correctly, but unfortunately they don’t give me the desired image.
I am rewriting the last code I wrote in Python here, maybe there is a way.
This code is written in a random and artificial way, but it still does not work, the real data has not plot a picture yet.
import numpy as np
import matplotlib.pyplot as plt
from sunpy.net import Fido, attrs
from astropy.time import Time
def get_hmi_data(start_time, end_time):
"""Get HMI data for the specified date."""
result = Fido.search(attrs.Time(start_time, end_time), attrs.Instrument('HMI'))
if result:
download = Fido.fetch(result)
return download
else:
print("No data found")
def generate_random_magnetic_field(num_lines=50, radius=1.0, noise_scale=0.1):
"""Generation of random magnetic field lines"""
lines = []
for _ in range(num_lines):
theta = np.random.uniform(0, 2 * np.pi)
phi = np.random.uniform(0, np.pi)
x = radius * np.sin(phi) * np.cos(theta)
y = radius * np.sin(phi) * np.sin(theta)
z = radius * np.cos(phi)
x += np.random.normal(0, noise_scale)
y += np.random.normal(0, noise_scale)
z += np.random.normal(0, noise_scale)
lines.append((x, y, z))
return lines
def plot_magnetic_field(lines):
"""plot magnetic field lines"""
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
for line in lines:
ax.scatter(line[0], line[1], line[2], c='b', marker='o')
ax.set_title('Artificial magnetic field lines')
ax.set_xlabel('X (Solar Radii)')
ax.set_ylabel('Y (Solar Radii)')
ax.set_zlabel('Z (Solar Radii)')
ax.set_xlim([-1.5, 1.5])
ax.set_ylim([-1.5, 1.5])
ax.set_zlim([-1.5, 1.5])
plt.show()
while True:
start_time_str = input("start date(formatYYYY-MM-DD): ")
end_time_str = input("end date (formatYYYY-MM-DD): ")
try:
start_time = Time(f'{start_time_str}T00:00:00')
end_time = Time(f'{end_time_str}T23:59:59')
break
except Exception as e:
print(f"Please try again {e}:Error in date format")
hmi_data = get_hmi_data(start_time, end_time)
if hmi_data:
num_lines = 100
radius = 1.0
lines = generate_random_magnetic_field(num_lines=num_lines, radius=radius)
plot_magnetic_field(lines)
else:
print("There is no data for this date")
در تاریخ چهارشنبه ۲۱ اوت ۲۰۲۴، ۱۹:۰۳ Nabil Freij via OpenAstronomy Community <notifications@openastronomy.discoursemail.com> نوشت:
Thanks for the tip, I’ve used the solved examples before and they worked correctly, but unfortunately they don’t give me the desired image.
What kind of image do you want to create?
This code is written in a random and artificial way, but it still does not work, the real data has not plot a picture yet.
I guess I would start by trying to piece together what you need and then when you check you have the previous step completed, add the next stage and continue.
My supervisor did not explain anything, but with my searches, I have obtained a series of images that should look like this, once the minimum state of the magnetic field and once the maximum state of the sun’s magnetic field.
در تاریخ چهارشنبه ۲۱ اوت ۲۰۲۴، ۲۰:۱۰ Nabil Freij via OpenAstronomy Community <notifications@openastronomy.discoursemail.com> نوشت:
So you want to replicate this image but from solar minimum and solar maximum?
I would strongly suggest asking your supervisor to explain more of the background.
I guess my other question is, do you need to create it yourself or if there are already images online that have this, will that be enough?
Yes , for solar maximum and solar minimum, for example in 2022 and 2023.
I’ll be honest, I want to do it myself. This is a part of my thesis and it is the final result of the efforts of all recent years.I have to make them myself, this is for my thesis result.
I have the code, but so far it has not given the output of the PFSS model.
Regarding the thesis supervisor, unfortunately I can’t do anything.
Is there an online map or figure of this model?
در تاریخ چهارشنبه ۲۱ اوت ۲۰۲۴، ۲۰:۳۰ Nabil Freij via OpenAstronomy Community <notifications@openastronomy.discoursemail.com> نوشت:
I have the code, but so far it has not given the output of the PFSS model.
This code being the one you gave above?
Do any of the examples I linked provide the correct image but the wrong date/data?
Is there an online map or figure of this model?
If you want a catalogue of those images, if you go to Sun In Time and go to the day you want, there might be a PFSS modelled image that you can use.
There is also this Interactive Solar Fieldline Viewer
I found several codes on GitHub as well as solved examples, some of which worked but had different output that was not relevant to my thesis. This is actually the last code I wrote today.
Thank you so much for every line of guidance, that’s the least I can say. I will try the links as well.
در تاریخ چهارشنبه ۲۱ اوت ۲۰۲۴، ۲۱:۴۶ Nabil Freij via OpenAstronomy Community <notifications@openastronomy.discoursemail.com> نوشت:
I found several codes on GitHub as well as solved examples, some of which worked but had different output that was not relevant to my thesis.
Do you have links to these codes?
When you say different output, was it just the style of the image or was it something else?
I have the file of one of them which I thought had the desired output.
I want to explain in detail, taking into account the number of sunspots, I have to find the day of maximum activity and minimum activity and put it in the code to finally have the solar PFSS model in two modes using the HMI tool.
در تاریخ چهارشنبه ۲۱ اوت ۲۰۲۴، ۲۲:۲۴ Nabil Freij via OpenAstronomy Community <notifications@openastronomy.discoursemail.com> نوشت:
(Attachment insunpygallery.py is missing)
Unfortunately that file was not uploaded.
Yes, unfortunately, I will write and send.
در تاریخ پنجشنبه ۲۲ اوت ۲۰۲۴، ۰۱:۲۵ Nabil Freij via OpenAstronomy Community <notifications@openastronomy.discoursemail.com> نوشت:
import astropy.constants as const
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import sunpy.map
from astropy.coordinates import SkyCoord
import pfsspy
from pfsspy import coords, tracing
from pfsspy.sample_data import get_gong_map
gong_fname = get_gong_map()
gong_map = sunpy.map.Map(gong_fname)
nrho = 35
rss = 2.5
pfss_in = pfsspy.Input(gong_map, nrho, rss)
def set_axes_lims(ax):
ax.set_xlim(0, 360)
ax.set_ylim(0, 180)
m = pfss_in.map
fig = plt.figure()
ax = plt.subplot(projection=m)
m.plot()
plt.colorbar()
ax.set_title(‘Input field’)
set_axes_lims(ax)
pfss_out = pfsspy.pfss(pfss_in)
ss_br = pfss_out.source_surface_br
Create the figure and axes
fig = plt.figure()
ax = plt.subplot(projection=ss_br)
Plot the source surface map
ss_br.plot()
Plot the polarity inversion line
ax.plot_coord(pfss_out.source_surface_pils[0])
Plot formatting
plt.colorbar()
ax.set_title(‘Source surface magnetic field’)
set_axes_lims(ax)
Get the radial magnetic field at a given height
ridx = 15
br = pfss_out.bc[0][:, :, ridx]
Create a sunpy Map object using output WCS
br = sunpy.map.Map(br.T, pfss_out.source_surface_br.wcs)
Get the radial coordinate
r = np.exp(pfss_out.grid.rc[ridx])
Create the figure and axes
fig = plt.figure()
ax = plt.subplot(projection=br)
Plot the source surface map
br.plot(cmap=‘RdBu’)
Plot formatting
plt.colorbar()
ax.set_title(‘B_{r} ’ + f’at r={r:.2f}’ + ‘r_{\\odot}’)
set_axes_lims(ax)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
tracer = tracing.FortranTracer()
r = 1.2 * const.R_sun
lat = np.linspace(-np.pi / 2, np.pi / 2, 8, endpoint=False)
lon = np.linspace(0, 2 * np.pi, 8, 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=pfss_out.coordinate_frame)
field_lines = tracer.trace(seeds, pfss_out)
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')
plt.show()
# sphinx_gallery_thumbnail_number = 4
در تاریخ پنجشنبه ۲۲ اوت ۲۰۲۴، ۰۹:۱۱ Zahra Atashkar <zahra.atashkar96@gmail.com> نوشت:
So this example code works for you and you want to now start to modify it for your date correct?
I can say yes, I want two maximums and two minimums, according to the number of recorded sunspots, I have a series of dates where the highest number is the maximum state and the lowest is the minimum state of the activity of the sun’s magnetic field.
در تاریخ پنجشنبه ۲۲ اوت ۲۰۲۴، ۲۰:۲۱ Nabil Freij via OpenAstronomy Community <notifications@openastronomy.discoursemail.com> نوشت:
So if that code works for you.
I would suggest that you should start by replacing the GONG image, to a different date and see what happens. The code should still work, you would just need to use the new filepath of the image.
To get GONG images, you should be able to search for it with Fido and download it that way.
If not, you can try the Virtual Solar Observatory via their website.
Thank you very much for your advice, so I will try to do it this way, maybe I will get the desired result. Thank you very much for every second and every word to guide me.
در تاریخ پنجشنبه ۲۲ اوت ۲۰۲۴، ۲۲:۱۱ Nabil Freij via OpenAstronomy Community <notifications@openastronomy.discoursemail.com> نوشت:
Hello.
Finally it worked and I was able to plot the maximum and minimum for 10 consecutive years.
در تاریخ پنجشنبه ۲۲ اوت ۲۰۲۴، ۲۲:۲۲ Zahra Atashkar <zahra.atashkar96@gmail.com> نوشت: