CLUSTER_Python_Programme/Bachelorthesis/Mollweide.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Apr 17 10:05:13 2026

@author: vieweg
"""

from rapid_fgm import *
from datetime import datetime, timedelta
import matplotlib.pyplot as plt

#--------------USER INPUT----------------------
TIME_PERIOD = (datetime(2002,6,27,00), datetime(2002,6,29,0))

TIME_POINT = datetime(2002,6,28,19,15)

SC_ID = 'C3'

e_idx = 6       #ch 0, ...,7

PARTICLE_TYPE = "electron"  #'electron' oder "ion"


#-----------------------------------------------
instrument = "IES" if PARTICLE_TYPE == "electron" else "IIMS"

r = RAPID(sc_id=SC_ID, time_period=TIME_PERIOD, particle_type=PARTICLE_TYPE)

data_rapid = r.load_rapid()
r.get_rapid_var(r.data_rapid)
print('rapid time:', r.time_rapid[0], r.time_rapid[-1])




flow = FLOW(
    sc_id=SC_ID,
    time_period=TIME_PERIOD,
    particle_type=PARTICLE_TYPE
)

flow.get_flow_var(flow.data_flow)
flow.load_flow()


f = FGM(SC_ID, time_period=TIME_PERIOD)
data_fgm = f.load_fgm()
f.get_fgm_var(f.data_fgm)

p = PITCH_CSA(SC_ID, time_period= TIME_PERIOD)
data_pad = p.load_pad()
pad_data = p.get_pad_var(p.data_pad)
#print(p.data_pad.cdf_info())

w = WORKSPACE() 
sync = w.sync_all_by_resolution(r, f, flow)

t_sync = sync['time']
rapid_sync = sync['rapid']
fgm_sync = sync['fgm']
sc2gse = sync['sc2gse']

sc_data = w.spacecrafts(
    sc_id=SC_ID,
    t_sync=t_sync,
    rapid_sync=rapid_sync,
    fgm_sync=fgm_sync,
    sc2gse=sc2gse,
    time_period=None)
print("len SC time:", len(w.SC_data['C3']['time']))



#font in Times New Roman
plt.rcParams['font.family'] = 'Times New Roman'
plt.rcParams['mathtext.fontset'] = 'stix'




w.rapid_n_sc2gse(SC_ID)


pitch_angle = w.calc_pitch_angle('C3')




def channel_label(e_idx):
    channel_phys = e_idx + 1
    energy = r.energy_table[e_idx]
    return energy, channel_phys


#  COLORS 
color1 ='#0070C0' 
color2 = '#63A375'


#  PAD

IES_energy, channel_phys = channel_label(e_idx)

pitch, flux, flux_sd = w.plot_flux_of_pitch(SC_ID, TIME_POINT, e_idx)
flux = np.ma.masked_where(flux <= 0, flux)
plt.figure(figsize=(10, 9))
plt.errorbar(
    pitch,
    flux,
    yerr=flux_sd,
    linestyle='none',
    marker='o',
    color='gray',
    ecolor='gray',
    label=r'Flux $\pm$ SD', alpha = 0.7
)
# Reference xerror bar
ref_x = 15                      

if len(flux) > 0:
    ref_y = np.nanmax(flux) * 0.9
else:
    ref_y = 1 


plt.errorbar(
    ref_x,
    ref_y,
    xerr=10, 
    fmt='o',
    color='black',
    capsize=4,
    label='10° Sector Width'
)

plt.xticks([10, 30, 50, 70, 90, 110, 130, 150, 170])
plt.xlabel('Pitch Angle / °')
plt.ylabel(r'Flux / $(\mathrm{cm^{2}\,s\,sr\,keV})^{-1}$')
plt.ylim(0,4000)

plt.title(
    f'{instrument} Pitch Angle Distribution - {SC_ID} - '
    f'Ch{channel_phys} ({IES_energy:.1f} keV) - {TIME_POINT}'
)

plt.legend()
plt.grid(True)
#plt.show()

#  Histogram Mean Flux

IES_energy, channel_phys = channel_label(e_idx)

bins, histo, shist, nhist = w.histo_flux_pitch(SC_ID, TIME_POINT, e_idx)
#durch anzahl pro bin teilen gewichtet mit fluss


print(w.SC_data[SC_ID]['time'][0], w.SC_data[SC_ID]['time'][-1])


#plot mean flux with errorbar
#plt.figure()
plt.errorbar(
    bins,
    histo,
    yerr=shist,
    color=color1,
    ecolor=color2,
    linestyle='none',
    linewidth = 4,
    markersize = 10,
    marker='^',
    capsize=6,
    label='Mean Flux ± SD'
)

plt.legend(fontsize = 14)

plt.title(
    f"Pitch Angle Distribution - {instrument} - {SC_ID} - "
    f"Ch{channel_phys} ({IES_energy:.1f} keV) - {TIME_POINT}", fontsize = 16
)

plt.grid(True)
plt.xticks([10, 30, 50, 70, 90, 110, 130, 150, 170], fontsize = 14)
plt.yticks(fontsize = 14)
plt.ylim(0,4000)
plt.ylabel(r'Flux / $(\mathrm{cm^{2}\,s\,sr\,keV})^{-1}$', fontsize = 16)
plt.xlabel(r'Pitch Angle / °', fontsize = 16)

timestamp = TIME_POINT.strftime("%Y%m%d_%H%M%S")
filename = 'C:/Users/lview/OneDrive/Dokumente/Desktop/C3_RAPID_20020628/191500/Flux_Pitch/'+ f'{SC_ID}_{instrument}_Ch'+f'{channel_phys}_'+f'{timestamp}'+'_FP_v2.pdf'
print(filename)
#plt.savefig(filename)
plt.show()


idx_pad = w.time_to_index(pad_data['time_pad'], TIME_POINT)
idx_sc  = w.time_to_index(w.SC_data[SC_ID]['time'], TIME_POINT)

#CSA PAD

flux_all = p.data_pad['PAD_Electron_Dif_flux__C3_CP_RAP_PAD_E3DD'][:]
flux_sd_all = p.data_pad['PAD_Electron_Dif_flux_SD__C3_CP_RAP_PAD_E3DD'][:]

pitch_csa = p.data_pad['Dimension_Pitch__C3_CP_RAP_PAD_E3DD'][:]
time_pad = p.data_pad['Time_tags__C3_CP_RAP_PAD_E3DD'][:]

import cdflib
time_pad_dt = cdflib.cdfepoch.to_datetime(time_pad)

idx_pad = w.time_to_index(time_pad_dt, TIME_POINT)

flux_csa = flux_all[idx_pad, e_idx, :]
flux_sd_csa = flux_sd_all[idx_pad, e_idx, :]
flux_sd_csa = flux_sd_all[idx_pad, e_idx, :]

plt.figure(figsize=(10,9))


plt.errorbar(
    pitch_csa,
    flux_csa,
    yerr=flux_sd_csa,
    fmt='o',
    color=color1,
    ecolor='green',
    capsize=5,
    markersize = 10,
    linewidth = 4,
    label='Flux ± SD'
)

# xerror separat
plt.errorbar(
    pitch_csa,
    flux_csa,
    xerr=10,
    fmt='none',          
    ecolor='black',
    capsize=5,
    linewidth = 4,
    label = '10° Sector Width')
  

plt.xticks([10, 30, 50, 70, 90, 110, 130, 150, 170], fontsize = 14)
plt.yticks(fontsize = 14)
plt.xlabel('Pitch Angle / °', fontsize = 16)
plt.ylabel(r'Flux / $(\mathrm{cm^{2}\,s\,sr\,keV})^{-1}$', fontsize = 16)
plt.ylim(0,4000)

plt.title(
    f'CSA {instrument} Pitch Angle Distribution - {SC_ID} - '
    f'Ch{channel_phys} ({IES_energy:.1f} keV) - {TIME_POINT}', fontsize = 16
)

plt.grid(True)
plt.legend(fontsize = 14)
timestamp = TIME_POINT.strftime("%Y%m%d_%H%M%S")
filename = 'C:/Users/lview/OneDrive/Dokumente/Desktop/C3_RAPID_20020628/191500/Flux_Pitch/'+ f'{SC_ID}_{instrument}_Ch'+f'{channel_phys}_'+f'{timestamp}'+'_FP_CSA.pdf'
print(filename)
plt.savefig(filename)
plt.show()