In [445]:
plotly.offline.init_notebook_mode(connected=True)
fig=plot_EM_loci_curves(fluxes,trmatrix=trmatrix,plotter='plotly')
fig.update_traces(line=dict(dash='dash'),opacity=0.6)
fig.for_each_trace(lambda trace: trace.update(name=trace.__getattribute__('name')+' predicted'))
fig.add_trace(go.Scatter(x=[np.log10(obs_params['T_spec_MK'][0]*1e6)],y=[obs_params['EM_spec_cm-3'][0]],error_x=errorx,error_y=errory,name='NuSTAR spectral fit',mode='markers',marker=dict(color='black',symbol='cross',line_width=2)))
for i,e in enumerate(evec[:-1]):
    fig.add_trace(go.Scatter(x=nustar_logt,y=o8rate[i]/o8tresp[:,i],name="%s-%sKeV" % (e,evec[i+1])))
fig.update_layout(yaxis_range=[40,45],xaxis_range=[5.5,8], title='AIA and NuSTAR loci curves, orbit %s' % orbit)

DEMreg with NuSTAR added

In [115]:
nlogt, nrate,ntresp=gen_nustar_tresp(eng_tr=[1.5,2.5,3.5], plot_loci=False)

/Users/wheatley/Documents/Solar/NuStar/specfiles/nu80610208001A06_chu12_N_sr_orbit8_2039_2042.arf /Users/wheatley/Documents/Solar/NuStar/specfiles/nu80610208001A06_chu12_N_sr_orbit8_2039_2042.rmf
In [132]:
lowE,highE=nustar_dem_prep(bl,tr)
lowE,highE #average counts/s in 1 px over entire flare peak
Out[132]:
(0.00076407596, 7.862101e-05)
In [337]:
lowEs,highEs=nustar_dem_prep(bl,tr,how='sum')
lowEs,highEs #average counts/s in entire area, over entire flare peak
Out[337]:
(0.05501347, 0.0056607127)
In [133]:
nt,nf,trmatrix,tresp_logt=gen_tresp_matrix(plot=True)

Want the lower temperatures for AIA as well, but f_vth used to calculate NuSTAR response shouldn't go below 6. So keep the AIA tresp_logt vector but make nustar response = 0 (or nan if that's accepted) for everything below 6

In [390]:
new_ntresp=np.zeros((len(tresp_logt),2))
new_ntresp[-40:,:]=ntresp 
#replace 0 with something very small
new_ntresp[new_ntresp == 0.0] = 1e-50

get AIA data over the same timerange, only the brightest parts

In [244]:
mask=make_contour_mask(211,submap=[bl20,tr20],tag=None,contour=[90],plot=True, diff=True)

1
In [245]:
np.product(mask.shape)- np.sum(mask) #number of pixels being used
Out[245]:
13
In [221]:
ts=dt.strptime('20200912_203700','%Y%m%d_%H%M%S')
te=dt.strptime('20200912_204200','%Y%m%d_%H%M%S')
In [253]:
dfaia.head()
Out[253]:
index data filenames fluxes maps rotated_coords timestamps wavelength cutout_shape
0 0 [[0.8685157398490652, 0.2211252322513833, 0.27... None 2.352308 [[0.86851574 0.22112523 0.2763082 ... 1.80209... (<SkyCoord (Helioprojective: obstime=2020-09-1... 2020-09-12T20:37:11.12 94 4422
1 1 [[0.5119616110509695, 1.7627540674154252, 1.48... None 2.552295 [[ 0.51196161 1.76275407 1.48023003 ... 0.4... (<SkyCoord (Helioprojective: obstime=2020-09-1... 2020-09-12T20:37:23.12 94 4422
2 2 [[2.942155809496954, 1.0679959877459109, 1.057... None 2.420161 [[2.94215581 1.06799599 1.05759489 ... 0.16255... (<SkyCoord (Helioprojective: obstime=2020-09-1... 2020-09-12T20:37:35.12 94 4422
3 3 [[0.6708007114998653, 0.5661923265307237, 1.82... None 2.602615 [[0.67080071 0.56619233 1.82239493 ... 0.22729... (<SkyCoord (Helioprojective: obstime=2020-09-1... 2020-09-12T20:37:47.12 94 4422
4 4 [[1.0508354614227233, 0.647897458954857, 0.347... None 2.319822 [[1.05083546 0.64789746 0.34733734 ... 0.67007... (<SkyCoord (Helioprojective: obstime=2020-09-1... 2020-09-12T20:37:59.14 94 4422
In [256]:
#now average over whole time period to geh the datavec
aia_datavec=[dfaia.where(dfaia.wavelength == w).dropna(how='all').fluxes.mean() for w in [94,131,171,193,211,335]]
aia_datavec #average DN/s in 1 pix inside 90% contour
Out[256]:
[2.374564492714406,
 10.088692964800845,
 331.7476385364175,
 403.56286623290305,
 82.49981603278998,
 1.7989615092177764]
In [258]:
#put it all together - data
aia_datavec.extend([lowE,highE])
aia_datavec
Out[258]:
[2.374564492714406,
 10.088692964800845,
 331.7476385364175,
 403.56286623290305,
 82.49981603278998,
 1.7989615092177764,
 0.00076407596,
 7.862101e-05]
In [303]:
datavec=np.array(aia_datavec)
In [392]:
#put it all together - response functions
tr_total=np.zeros((101,8))
tr_total[:,:6]=trmatrix#_trim
tr_total[:,6:]=new_ntresp
In [416]:
#run the DEM
temps=np.logspace(5.7,7.5,num=40)  #or use Mark's bins (21?) og 35
dem_norm0=dem_norm_guess(temps,1,1,40)
edata=generate_errors(0,0,6,aia_datavec[:-2]) #what are the nustar errors? do sqrt for now... but these are larger numbers than the data
edata_ns=np.sqrt([lowE,highE])
In [417]:
edata_total=np.zeros(8)
edata_total[:6]=edata
edata_total[6:]=edata_ns
edata_total
Out[417]:
array([2.48604537e+00, 4.16137033e+00, 3.83970689e+01, 4.49245320e+01,
       1.19195923e+01, 1.58229476e+00, 2.76419241e-02, 8.86684842e-03])
In [418]:
#using summed data
edatas=generate_errors(0,0,6,aia_datavecs[:-2]) #what are the nustar errors? do sqrt for now... but these are larger numbers than the data
edata_nss=np.sqrt([lowEs,highEs])
edata_totals=np.zeros(8)
edata_totals[:6]=edatas
edata_totals[6:]=edata_nss
edata_totals
Out[418]:
array([8.48079194e+00, 1.91489617e+01, 4.36849252e+02, 5.29420266e+02,
       1.11568283e+02, 4.51423411e+00, 2.34549507e-01, 7.52377063e-02])
In [419]:
dem,edem,elogt,chisq,dn_reg=dn2dem_pos(datavec,edata_total,tr_total,tresp_logt,temps,dem_norm0=dem_norm0)

Executing in serial
In [420]:
#using sum
dems,edems,elogts,chisqs,dn_regs=dn2dem_pos(datavecs,edata_totals,tr_total,tresp_logt,temps,dem_norm0=dem_norm0)

Executing in serial
In [421]:
#using sum, without nustar
dems0,edems0,elogts0,chisqs0,dn_regs0=dn2dem_pos(datavecs[:6],edata_totals[:6],trmatrix,tresp_logt,temps,dem_norm0=dem_norm0)

Executing in serial
In [422]:
#without nustar
dem0,edem0,elogt0,chisq0,dn_reg0=dn2dem_pos(datavec[:6],edata_total[:6],trmatrix,tresp_logt,temps,dem_norm0=dem_norm0)

Executing in serial
In [433]:
fig=go.Figure()
fig.add_trace(go.Scatter(x=np.log10(temps),y=dem,error_x=dict(type='data',array=elogt),error_y=dict(type='data',array=edem),name='DEM with NuSTAR (mean)'))
fig.add_trace(go.Scatter(x=np.log10(temps),y=dem0,error_x=dict(type='data',array=elogt0),error_y=dict(type='data',array=edem0),name='DEM without NuSTAR (mean)'))
fig.add_trace(go.Scatter(x=np.log10(temps),y=dems,error_x=dict(type='data',array=elogts),error_y=dict(type='data',array=edems),name='DEM with NuSTAR (sum)'))
fig.add_trace(go.Scatter(x=np.log10(temps),y=dems0,error_x=dict(type='data',array=elogts0),error_y=dict(type='data',array=edems0),name='DEM without NuSTAR (sum)'))
fig.update_layout(yaxis_type='log',xaxis_title='log T',yaxis_title='DEM (cm<sup>5</sup> K<sup>-1</sup>)', title="DEMreg result")
fig.update_layout(yaxis = dict(showexponent = 'all',exponentformat = 'e'))
In [429]:
xvec=['094','131','171','193','211','335','1.5-2.5','2.5-3.5']
fig=go.Figure()
fig.add_trace(go.Scatter(x=xvec,y=dn_reg-datavec,name='with NuSTAR (mean)'))
fig.add_trace(go.Scatter(x=xvec[:-2],y=dn_reg0-datavec[:-2],name='without NuSTAR (mean)'))
fig.add_trace(go.Scatter(x=xvec,y=dn_regs-datavecs,name= 'with NuSTAR (sum)'))
fig.add_trace(go.Scatter(x=xvec[:-2],y=dn_regs0-datavecs[:-2],name='without NuSTAR (sum)'))
fig.update_layout(xaxis_title='wavelength',yaxis_title='Difference',title='dn_reg - data')
#fig.update_layout(yaxis = dict(showexponent = 'all',exponentformat = 'e'))