In [ ]:
%matplotlib inline
In [112]:
# imports
from desispec import bootcalib as desiboot
from astropy.io import fits
from astropy.stats import sigma_clip
import numpy as np
from astropy.modeling import models, fitting
from xastropy.xutils import afits as xafits
In [2]:
# Read flat
flat = fits.open('/Users/xavier/DESI/Wavelengths/pix-b0-00000001.fits')[0].data
In [3]:
reload(desiboot)
xpk, ypos, cut = desiboot.find_fiber_peaks(flat)#,debug=True)
In [4]:
# Check
plt.clf()
plt.figure(figsize=(16,7))
xplt = np.arange(cut.size)
plt.plot(xplt,cut, 'k-')
plt.plot(xpk, cut[xpk],'go')
plt.xlim(100,500)
plt.show()
In [5]:
reload(desiboot)
xset, xerr = desiboot.trace_crude_init(flat,xpk[0:50],ypos)
In [6]:
yval = np.arange(4096)
ii=4
xval = xset[:,ii]
gdval = xerr[:,ii] < 999.
dfit0 = xafits.func_fit(yval[gdval],xval[gdval],'legendre',6)
In [7]:
fitv = xafits.func_val(yval,dfit0)
In [8]:
# Fit
xdb.xplot(yval,fitv, xtwo=yval,ytwo=xset[:,ii],mtwo='+')
In [ ]:
# Residuals
xdb.xplot(yval[gdval],fitv[gdval]-xval[gdval], scatter=True)
In [9]:
reload(desiboot)
xfit, fdicts = desiboot.fit_traces(xset,xerr)#[:,0:5],xerr[:,0:5])
In [10]:
reload(desiboot)
desiboot.fiber_trace_qa(flat,xfit)
In [52]:
reload(desiboot)
gauss = desiboot.fiber_gauss(flat,xfit,xerr)#,verbose=True)#,debug=True)
In [25]:
fiber = np.arange(gauss.size)
gfdict,mask = xafits.iter_fit(fiber, gauss, 'polynomial', 2)
In [26]:
gfdict
Out[26]:
In [27]:
plt.clf()
plt.scatter(fiber,gauss)
plt.plot(fiber, xafits.func_val(fiber,gfdict),'k-')
#plt.xlim(100,500)
plt.show()
In [48]:
# Load
arc = fits.open('/Users/xavier/DESI/Wavelengths/pix-b0-00000000.fits')[0].data
In [67]:
sv_aspec = []
In [340]:
ii=0
In [341]:
xpix_img = np.outer(np.ones(arc.shape[0]),np.arange(arc.shape[1]))
mask = np.zeros_like(arc,dtype=int)
iy = np.arange(arc.shape[0],dtype=int)
box_radius = 2
ixt = np.round(xfit[:,ii]).astype(int)
for jj,ibox in enumerate(range(-box_radius,box_radius+1)):
ix = ixt + ibox
mask[iy,ix] = 1
In [342]:
dx_img = xpix_img - np.outer(xfit[:,ii],np.ones(arc.shape[1]))
g_init = models.Gaussian1D(amplitude=1., mean=0., stddev=gauss[ii])
psf = mask * g_init(dx_img)
#xdb.ximshow(psf)
In [343]:
arc_spec = np.sum(psf*arc,axis=1) / np.sum(psf,axis=1)
In [84]:
sv_aspec.append(arc_spec)
In [85]:
#xdb.xplot(np.arange(arc_spec.size), sv_aspec[0], sv_aspec[1], sv_aspec[2])
In [29]:
# Dispersion
import desimodel.io
In [87]:
desi_psf = desimodel.io.load_psf('b')
In [88]:
desi_psf.wdisp(0,4500.) # fiber, wavelength
Out[88]:
In [619]:
wave0 = desi_psf.wavelength(0,np.arange(desi_psf.npix_y))
In [135]:
np.max(wave0)
Out[135]:
In [132]:
med_bdisp = np.median(desi_psf.wdisp(0,wave0))
med_bdisp
Out[132]:
In [626]:
np.median(np.abs(wave0-np.roll(wave0,1)))
Out[626]:
In [623]:
xdb.xplot(wave0,desi_psf.wdisp(0,wave0))
In [231]:
desi_psf.npix_x, desi_psf.npix_y
Out[231]:
In [624]:
poly_fit = xafits.func_fit(wave0,np.arange(desi_psf.npix_y), 'polynomial',2,xmin=0.,xmax=1.)
poly_fit
Out[624]:
In [625]:
xdb.xplot(wave0, np.arange(desi_psf.npix_y),xafits.func_val(wave0,poly_fit))
In [348]:
# Set flux threshold
cut = arc_spec
cut_mask = sigma_clip(cut,sig=4.)
rms = np.std(cut_mask)
thresh = 10*rms
print("thresh = {:g}".format(thresh))
gdp = cut > thresh
print("Pixels exceeding = {:d}".format(np.sum(gdp)))
In [349]:
# Roll to find peaks (simple algorithm)
nwidth = 5
nstep = nwidth // 2
for kk in xrange(-nstep,nstep):
if kk < 0:
test = np.roll(cut,kk) < np.roll(cut,kk+1)
else:
test = np.roll(cut,kk) > np.roll(cut,kk+1)
# Compare
gdp = gdp & test
In [350]:
xpk = np.where(gdp)[0]
In [351]:
plt.clf()
yspec = np.log10(np.maximum(cut,1))
xplt = np.arange(cut.size)
plt.plot(xplt,yspec,'b-')
plt.scatter(xpk,yspec[xpk],color='green')
plt.ylim(0.,np.max(yspec)*1.05)
plt.xlim(0,cut.size)
plt.show()
In [127]:
sys.path.append(os.path.abspath(os.getenv('PYPIT')+"/src"))
import ararclines as pypit_alines
import ararc as pypit_arc
In [627]:
arcparam = dict(llist=None,
disp=0.592, # Ang (unbinned)
b1=0., # Pixel fit term (binning independent)
b2=0., #1.09530252e-05/arc.shape[0], # Pixel fit term
wvmnx=[3500.,6000.],# Pixel fit term
disp_toler=0.1, # 10% tolerance
match_toler=5., # Matcing tolerance (pixels)
func='legendre', # Function for fitting
n_first=1, # Order of polynomial for first fit
n_final=4, # Order of polynomial for final fit
nsig_rej=2., # Number of sigma for rejection
nsig_rej_final=3.0, # Number of sigma for rejection (final fit)
Nstrong=15)
arcparam['b1']= 1./arcparam['disp']/arc.shape[0]
aparm = arcparam
aparm
Out[627]:
In [796]:
reload(pypit_alines)
#llist = pypit_alines.load_arcline_list(None,None,['CdI','ArI','NeI','HgI'],wvmnx=aparm['wvmnx'])
llist = pypit_alines.load_arcline_list(None,None,['CdI','ArI','HgI','NeI'],wvmnx=aparm['wvmnx'])
In [629]:
arcparam['llist'] = llist
In [797]:
llist['wave'][llist['Ion']=='NeI']
Out[797]:
In [798]:
tampl = arc_spec[xpk]
tcent = xpk
# Need to improve centroid to sub-pixel (flux-weight)
In [799]:
# Generate dpix pairs
nlist = len(llist)
dpix_list = np.zeros((nlist,nlist))
for kk,row in enumerate(llist):
dpix_list[kk,:] = arc.shape[0]*(aparm['b1']*(np.array(row['wave'] - llist['wave'])) + aparm['b2']*np.array(row['wave']**2 - llist['wave']**2) )
In [800]:
dpix_list[8,:]
Out[800]:
In [801]:
dpix_list[13,:]
Out[801]:
In [802]:
# Pixel pairs for the strongest N lines
srt = np.argsort(tampl)
idx_str = srt[-aparm['Nstrong']:]
idx_str.sort()
dpix_obs = np.zeros((aparm['Nstrong'],aparm['Nstrong']))
for kk,idx in enumerate(idx_str):
dpix_obs[kk,:] = np.array(tcent[idx] - tcent[idx_str]) # pixels
In [803]:
dpix_obs[8,:]
Out[803]:
In [804]:
# Match up (ugly loops)
ids = np.zeros(aparm['Nstrong'])
ids2 = np.zeros(aparm['Nstrong'])
ions = np.array(['12345']*aparm['Nstrong'])
for kk in range(aparm['Nstrong']):
med_off = np.zeros(nlist)
mean_off = np.zeros(nlist)
for ss in range(nlist):
dpix = dpix_list[ss]
min_off = []
for jj in range(aparm['Nstrong']):
min_off.append(np.min(np.abs(dpix_obs[kk,jj]-dpix)))
min_off.sort()
mean_off[ss] = np.mean(min_off[:-1])
med_off[ss] = np.median(min_off)
#xdb.set_trace()
# if (kk == 1) & ((ss == 13) or (ss==8)):
# xdb.set_trace()
# Set by minimum
idm = np.argmin(med_off)
idm2 = np.argmin(mean_off)
print(kk,mean_off[idm2])
ids2[kk] = llist['wave'][idm]
ids[kk] = llist['wave'][idm2]
ions[kk] = llist['Ion'][idm]
In [805]:
xdb.xpcol(ids,ids2)
In [806]:
# Calculate disp of the strong lines
disp_str = np.zeros(aparm['Nstrong'])
for kk in range(aparm['Nstrong']):
disp_str[kk] = np.median( (ids[kk]-ids)/
(tcent[idx_str[kk]]-tcent[idx_str]) )
# Consider calculating the RMS with clipping
gd_str = np.where( np.abs(disp_str-aparm['disp'])/aparm['disp'] < aparm['disp_toler'])[0]
print('Found {:d} lines within the disperesion threshold'.format(len(gd_str)))
if len(gd_str) < 5:
raise ValueError('Insufficient lines to auto-fit.')
In [807]:
xdb.xpcol(np.arange(len(ids)),ids,disp_str)
In [808]:
#Plot
plt.clf()
xplt = np.arange(arc_spec.size)
plt.plot(xplt,arc_spec)
# Label
for gds in gd_str:
plt.text(tcent[idx_str[gds]], tampl[idx_str[gds]], '{:s} {:g}'.format(ions[gds],ids[gds]), rotation=90., va='bottom')
plt.show()
In [415]:
# Fit
#Plot
plt.clf()
plt.scatter(ids[gd_str],tcent[idx_str[gd_str]],color='blue')
plt.plot(wave0,np.arange(wave0.size),'g-')
plt.show()
In [642]:
# Residual
resid = tcent[idx_str[gd_str]] - xafits.func_val(ids[gd_str], poly_fit)
plt.clf()
plt.scatter(ids[gd_str],resid,color='blue')
plt.show()
In [412]:
ids[gd_str]
Out[412]:
In [765]:
llist.sort('wave')
nline = len(llist)
In [583]:
llist['wave']
Out[583]:
In [ ]:
# Generate
#3610, 3650., 3654
line_triplets = []
for kk in range(1,nline-1):
#print(kk)
# Setup array
triplets = np.zeros( (3,kk,nline-kk-1))
wv2 = llist['wave'][kk]
# Loop on LHS
for jj in range(kk):
wv1 = llist['wave'][jj]
wvmm = llist['wave'][kk+1:]
# if kk == 2:
# xdb.set_trace()
triplets[0,jj,:] = (wv2-wv1)/(wvmm-wv1)
triplets[1,jj,:] = (wv2-wv1)/(wvmm-wv2)
triplets[2,jj,:] = (wvmm-wv2)/(wvmm-wv1)
# if kk == 34:
# if np.sum(np.abs(triplets[1,jj,:]-8)<1e-5)>0:
# xdb.set_trace()
#triplets[3,jj,:] = wvmm-wv1 # dlamb
# Enforce 0-1 range
#if kk == 1:
# xdb.set_trace()
#inverse = triplets > 1.
#triplets[inverse] = 1./triplets[inverse]
# Reshape
# if kk == 34:
# xdb.set_trace()
triplets = np.reshape(triplets,(3,triplets.shape[1]*triplets.shape[2]))
# Sort
#triplets = np.sort(triplets, axis=0)
# Append
#xdb.set_trace()
line_triplets.append(triplets)
In [484]:
tcent[0:3]
Out[484]:
In [586]:
ii = 1
triplet_obs = [ float(tcent[ii]-tcent[ii-1])/float(tcent[ii+1]-tcent[ii-1]),
float(tcent[ii]-tcent[ii-1])/float(tcent[ii+1]-tcent[ii]),
float(tcent[ii+1]-tcent[ii])/float(tcent[ii+1]-tcent[ii-1])]
print(tcent[ii-1:ii+2])
triplet_o = np.array(triplet_obs)
# inverse
#inverse = triplet_o > 1.
#triplet_o[inverse] = 1./triplet_o[inverse]
#
#triplet_obs = np.sort(triplet_o)
triplet_o
Out[586]:
In [ ]:
# Loop on lines
all_diff = []
for kk,ltrip in enumerate(line_triplets):
diff = np.abs(ltrip.T - triplet_o)
min_diff = np.min(np.sum(diff,axis=1))
imin = np.argmin(np.sum(diff,axis=1))
if min_diff < 0.03:
pass
#dlamb = 1./(ltrip[3,imin]*(tcent[ii+1]-tcent[ii-1]))
#print('wave={:g}, min={:g}, dlamb={:g}'.format(llist['wave'][kk+1],min_diff))#,dlamb))
#all_diff.append(ltrip)
#xdb.set_trace()
#print(kk,min_diff)
#if kk == 33:
# xdb.set_trace()
In [570]:
ltrip = line_triplets[0]
In [558]:
diff = np.abs(ltrip[:,0:3] - triplet_o)
In [559]:
np.min(np.sum(diff,axis=1))
Out[559]:
In [552]:
xdb.xplot(arc_spec)
In [809]:
llist.sort('wave')
nline = len(llist)
nline
Out[809]:
In [810]:
npair = np.sum(np.arange(nline))
npair
Out[810]:
In [811]:
Dlamb = np.zeros((3,npair)) # index, index, dlamb
count = 0
for kk in range(nline-1):
jidx = np.arange(kk+1,nline)
cidx = np.arange(jidx.size)
dl = llist['wave'][jidx] - llist['wave'][kk]
# Fill
Dlamb[0,count+cidx] = kk
Dlamb[1,count+cidx] = jidx
Dlamb[2,count+cidx] = dl
count += cidx.size
In [812]:
toler = 0.2
In [813]:
dlamb = (wave0[-1]-wave0[0])/wave0.size
dlamb
Out[813]:
In [814]:
#
Dobs = (tcent[-1]-tcent[0])*dlamb
Dobs
Out[814]:
In [815]:
#Normalized pattern of tcent
xtcent = (tcent[1:-1]-tcent[0])/float(tcent[-1]-tcent[0])
xtcent
Out[815]:
In [816]:
# All matches within tolerance
gdp = np.where(np.abs(Dobs-Dlamb[2,:])/Dobs < toler)[0]
gdp.size
Out[816]:
In [817]:
#
mtch_wv = np.zeros( (tcent.size, gdp.size))
mtch_idx = np.zeros( (tcent.size, gdp.size),dtype=int)
In [818]:
# Loop on gdp
metric = np.zeros(gdp.size)
for qq,igd in enumerate(gdp):
# Endpoints
i0,i1 = Dlamb[0,igd], Dlamb[1,igd]
mtch_wv[0,qq] = llist['wave'][i0]
mtch_wv[-1,qq] = llist['wave'][i1]
mtch_idx[0,qq] = i0
mtch_idx[-1,qq] = i1
# Normalized pattern of lines
xline = (llist['wave'][i0+1:i1]-llist['wave'][i0])/(llist['wave'][i1]-llist['wave'][i0])
# Minimize tcent
for jj,xt in enumerate(xtcent):
diff = np.abs(xt-xline)
imin = np.argmin(diff)
# Increment and save
metric[qq] += diff[imin]
mtch_wv[jj+1,qq] = llist['wave'][i0+imin+1]
mtch_idx[jj+1,qq] = i0+imin+1
In [819]:
xdb.xhist(metric)
In [777]:
best_met = np.argmin(metric)
best_met
Out[777]:
In [820]:
mtch_wv[:,best_met]
Out[820]:
In [821]:
mtch_idx[:,best_met]
Out[821]:
In [822]:
srt = np.argsort(metric)
In [823]:
srt[0:5]
Out[823]:
In [824]:
gdwv = np.median(mtch_wv[:,srt[0:5]],axis=1)
gdwv
Out[824]:
In [825]:
gdidx = np.median(mtch_idx[:,srt[0:5]],axis=1)
gdidx
Out[825]:
In [846]:
#Plot with labels
plt.clf()
xplt = np.arange(arc_spec.size)
plt.plot(xplt,arc_spec)
# Label
for ii in range(tcent.size):
plt.text(tcent[ii], tampl[ii], '{:g}'.format(gdwv[ii]), rotation=90., va='bottom')
plt.show()
In [845]:
xdb.xplot(wave0,arc_spec)
In [840]:
iq=14 # Marked line
In [841]:
quint = tcent[iq+np.arange(5)-2]
quint
Out[841]:
In [842]:
guesses = np.arange(gdidx[iq-3],gdidx[iq+1]+1)
qmetric = np.ones(guesses.size) * 999.
#
Dlm1 = dlamb*(tcent[iq]-tcent[iq-1])
Dlp1 = dlamb*(tcent[iq+1]-tcent[iq])
#
x1 = (tcent[iq-1]-tcent[iq-2])/float(tcent[iq]-tcent[iq-2])
x2 = (tcent[iq]-tcent[iq-1])/float(tcent[iq+1]-tcent[iq-1])
x3 = (tcent[iq+1]-tcent[iq])/float(tcent[iq+2]-tcent[iq])
for ss,icen in enumerate(guesses):
wcen = llist['wave'][icen]
print('wcen = {:g}'.format(wcen))
# Find possible i-1
Dlm = wcen-llist['wave'][0:icen]
gdm = np.where(np.abs(Dlm-Dlm1)/Dlm1<toler)[0]
# Find possible i+1
Dlp = llist['wave'][icen+1:]-wcen
gdp = np.where(np.abs(Dlp-Dlp1)/Dlp1<toler)[0]
#
nposs = len(gdm)*len(gdp)
if (nposs == 0):
continue
# Ugly loops
sub_metric = np.zeros(nposs)
count=0
for igdm in gdm:
# 012
lx1 = (llist['wave'][igdm]-llist['wave'][0:igdm])/(wcen-llist['wave'][0:igdm])
sub_metric[count+np.arange(gdp.size)] += np.min(np.abs(lx1-x1))
for igdp in gdp:
# 123
sub_metric[count] += np.abs(x2-(wcen-llist['wave'][igdm])/(llist['wave'][icen+1+igdp]-llist['wave'][igdm]))
count += 1
count=0
for igdp in gdp:
# 234
lx3 = (llist['wave'][icen+1+igdp]-wcen)/(llist['wave'][icen+2+igdp:]-wcen)
sub_metric[count*gdm.size+np.arange(gdm.size)] += np.min(np.abs(lx3-x3))
#xdb.set_trace()
qmetric[ss] = np.min(sub_metric)
In [843]:
xdb.xpcol(llist['wave'][guesses.astype(int)], qmetric)
In [844]:
xdb.xplot(arc_spec)
In [961]:
from astropy.modeling import models, fitting
import copy
In [893]:
# For fitting
fitter = fitting.LevMarLSQFitter()
In [891]:
# Setup model
poly = models.Polynomial1D(2)
poly
Out[891]:
In [895]:
# Set bounds
poly.c1.bounds = (6.38522295e-01*0.8, 6.38522295e-01*1.2)
In [870]:
Out[870]:
In [872]:
xdb.xplot(wave0,arc_spec)
In [874]:
HgI = [4046.57, 4077.84, 4358.34, 5460.75, 5769.598]
CdI = [3610.51, 3650.157, 4678.15, 4799.91, 5085.82]
NeI = [5085.822, 5881.895, 5944.834]
gd_lines = np.array(HgI + CdI + NeI)
gd_lines.sort()
gd_lines
Out[874]:
In [878]:
dlamb, toler
Out[878]:
In [929]:
rms_dicts = []
## Assign a line to focus on
#wmark = 4358.34
wmark = 4799.91
#wmark = 5085.822
icen = np.where(np.abs(gd_lines-wmark)<1e-3)[0]
icen = icen[0]
##
guesses = 8 + np.arange(-4,7)
for guess in guesses:
# Setup
tc = float(tcent[guess])
print('tc = {:g}'.format(tc))
rms_dict = dict(tc=tc,guess=guess)
## Match to i-2 line
line_m2 = gd_lines[icen-2]
Dm2 = wmark-line_m2
mtm2 = np.where(np.abs((tc-tcent)*dlamb - Dm2)/Dm2 < toler)[0]
if len(mtm2) == 0:
print('No im2 lines found for guess={:g}'.format(tc))
continue
# Match to i+2 line
line_p2 = gd_lines[icen+2]
Dp2 = line_p2-wmark
mtp2 = np.where(np.abs((tcent-tc)*dlamb - Dp2)/Dp2 < toler)[0]
if len(mtp2) == 0:
print('No ip2 lines found for guess={:g}'.format(tc))
continue
#
all_guess_rms = [] # One per set of guess, of course
# Loop on i-2 line
for imtm2 in mtm2:
#
if imtm2==(icen-1):
print('No im1 lines found for guess={:g}'.format(tc))
continue
# Setup
tcm2 = float(tcent[imtm2])
xm1_wv = (gd_lines[icen-1]-gd_lines[icen-2])/(wmark-gd_lines[icen-2])
# Best match
xm1 = (tcent-tcm2)/(tc-tcm2)
itm1 = np.argmin(np.abs(xm1-xm1_wv))
# Now loop on ip2
for imtp2 in mtp2:
guess_rms = dict(guess=guess, im1=itm1,im2=imtm2, rms=999., ip1=None, ip2=imtp2)
all_guess_rms.append(guess_rms)
if imtp2==(icen-1):
print('No ip1 lines found for guess={:g}'.format(tc))
continue
#
tcp2 = float(tcent[imtp2])
xp1_wv = (gd_lines[icen+2]-gd_lines[icen+1])/(gd_lines[icen+2]-wmark)
# Best match
xp1 = (tcp2-tcent)/(tcp2-tc)
itp1 = np.argmin(np.abs(xp1-xp1_wv))
guess_rms['ip1'] = itp1
# Time to fit
xval = np.array([tcm2,float(tcent[itm1]),tc,float(tcent[itp1]),tcp2])
wvval = gd_lines[icen-2:icen+3]
pfit = xafits.func_fit(wvval,xval,'polynomial',2,xmin=0.,xmax=1.)
#parm = fitter(poly,wvval,xval)
# Clip one here and refit
# RMS (in pixel space)
rms = np.sqrt(np.sum((xval-xafits.func_val(wvval,pfit))**2)/xval.size)
guess_rms['rms'] = rms
# Save fit too
guess_rms['fit'] = pfit
#print('rms = {:g}'.format(rms))
#if np.abs(rms-4.93) < 1e-1:
# xdb.set_trace()
# Take best RMS
if len(all_guess_rms) > 0:
all_rms = np.array([idict['rms'] for idict in all_guess_rms])
imn = np.argmin(all_rms)
rms_dicts.append(all_guess_rms[imn])
#xdb.set_trace()
In [937]:
all_rms = np.array([idict['rms'] for idict in rms_dicts])
imin = np.argmin(all_rms)
best_dict = rms_dicts[imin]
best_dict
Out[937]:
In [963]:
wvval0 = list(gd_lines[icen-2:icen+3])
wvval0
Out[963]:
In [964]:
xval0 = []
for key in ['im2','im1','guess','ip1','ip2']:
xval0.append(float(tcent[best_dict[key]]))
xval0
Out[964]:
In [971]:
ilow = icen-2
ihi = icen+2
print(ilow, ihi)
#
npoly=2
pos=True
wvval = copy.deepcopy(wvval0)
xval = copy.deepcopy(xval0)
while ((ilow > 0) or (ihi < gd_lines.size-1)):
# index to add (step on each side)
if pos:
ihi += 1
inew = ihi
pos=False
else:
ilow -= 1
inew = ilow
pos=True
if ilow < 0:
continue
if ihi > (gd_lines.size-1):
continue
#print('inew={:d}'.format(inew))
# New line
new_wv = gd_lines[inew]
wvval.append(new_wv)
wvval.sort()
# newx
newx = xafits.func_val(new_wv,best_dict['fit'])
# Match and add
imin = np.argmin(np.abs(tcent-newx))
newtc = float(tcent[imin])
xval.append(newtc)
xval.sort()
#print('Added x={:g} at wv={:g}'.format(newtc,new_wv))
#xdb.set_trace()
# Fit (reject 1)
if len(xval) > 7:
npoly = 3
new_fit = xafits.func_fit(np.array(wvval),np.array(xval),'polynomial',npoly,xmin=0.,xmax=1.)
best_dict['fit'] = new_fit
# RMS
resid2 = (np.array(xval)-xafits.func_val(np.array(wvval),best_dict['fit']))**2
rms = np.sqrt(np.sum(resid2)/len(xval))
print('rms = {:g}'.format(rms))
#print(resid2)
#print(xval)
#print(wvval)
In [972]:
best_dict
Out[972]:
In [967]:
xdb.xplot(arc_spec)
In [ ]: