; ---------------------------------------------------------------------- ; Extract sky spectrum and compute skyline shift in pixels. ; Used in nsx_extract_spectrum() routine. ; IMGSKY : Sky data image (A image for both A only and A-B data). ; Loads .spsky[], .spwav[], and .spdsp[]. The 'spwav' values are shifted ; by a median pixel value (written to log file). ; ;> Fills SPX.spwav, spdsp, spsky ;> writes -skyres.tbl pro nsx_compute_skyline_shift, IMGSKY, SPX, nsxdir, nsxout common nsx ssnso = intarr(3000) skynum = intarr(9) xx = dblarr(3000) yy = dblarr(3000) arr = dblarr(1000) skywave = dblarr(9,100) skywgt = dblarr(9,100) sscent = dblarr(1000) sscol = dblarr(1000) ssresi = dblarr(1000) sspeak = dblarr(1000) ssfwhm = dblarr(1000) sswave = dblarr(1000) xcc = dblarr(3000) ycc = dblarr(3000) skytp = dblarr(3000) xx8 = dblarr(100) yy8 = dblarr(100) ww8 = dblarr(100) coef8 = dblarr(9) print, "Compute skyline shift." ; For each order, load wavelength scale and dispersion. ;> sets SPX.spwav and SPX.spdsp for nso=3,7 do begin if (nso eq 7) then ecol=IMGSKY.nc/2 else ecol=IMGSKY.nc for icol=0,ecol-1 do begin xv = double(icol) - WSC[nso].xoff SPX[nso].spwav[icol] = cpolyval(WSC[nso].order+1,WSC[nso].coef,xv) endfor ; Compute dispersion (angstroms per pixel) (about -2.8 ang/pix). for icol=0,ecol-1 do begin if (icol eq 0 ) then begin SPX[nso].spdsp[icol] = SPX[nso].spwav[1] - SPX[nso].spwav[0] endif else begin if (icol eq ecol-1) then begin SPX[nso].spdsp[icol] = SPX[nso].spwav[ecol-1] - SPX[nso].spwav[ecol-2] endif else begin SPX[nso].spdsp[icol] =(SPX[nso].spwav[icol+1] - SPX[nso].spwav[icol-1])/2. endelse endelse endfor ; icol endfor ; nso ; Extract sky spectrum from corrected background image for each order (in object window). ;> sets SPX.spsky for nso=3,7 do begin if (nso eq 7) then ecol=IMGSKY.nc/2 else ecol=IMGSKY.nc jjoff = ((7-nso) * 200) rblo = double(jjoff) rbhi = double(jjoff) + SPX[nso].numpro for icol=0L,ecol-1 do begin rb1 = double(jjoff) + (SPX[0].asp1[0] / ARCSEC_PER_PIXEL) if (rb1 lt rblo) then rb1 = rblo if (rb1 gt rbhi) then rb1 = rbhi rb2 = double(jjoff) + (SPX[0].asp2[0] / ARCSEC_PER_PIXEL) if (rb2 lt rblo) then rb2 = rblo if (rb2 gt rbhi) then rb2 = rbhi sum = nsx_fractional_pixel_rb( IMGSKY.nc, IMGSKY.bckimg, rb1, rb2, icol ) npx = ABS((rb2 - rb1)) if (npx gt 0.) then SPX[nso].spsky[icol] = sum/npx else SPX[nso].spsky[icol] = 0. SPX[nso].spsky[icol] = SPX[nso].spsky[icol] / IMGSKY.exptime if nso eq 0 then stop ;> DEBUG endfor ;> icol endfor ;> nso ; Read in sky line wavelength data. for nso=3,7 do begin wrd = string(format='(A,"cal/skylines_nsx",I1,".tbl")',nsxdir,nso) ;> use readcol readcol, wrd, _skywave, _skywgt, count=kk, skip=1, for='D,D', /silent skywave[nso,0:kk-1] = _skywave skywgt[nso,0:kk-1] = _skywgt skynum[nso]=kk endfor ; Cross correlation .. narr=0 for nso=3,7 do begin if (nso eq 7) then ecol=IMGSKY.nc/2 else ecol=IMGSKY.nc for icol=0,3000-1 do skytp[icol]=0. wrd = string(format='(A,"cal/sky_template",I1,".tbl")',nsxdir,nso) readcol, wrd, _icol, _skytp, count=count, skip=1, for='L,D,D', /silent skytp[_icol] = _skytp print, 'skytp[1024] ', skytp[1024] ;> DEBUG ; sprintf(wrd,"%s/cal/sky_template%d.tbl",nsxdir,nso) ; infu = fopen_read(wrd) ; while (fgetline(line,infu)) { if (line[0] != '|') { ; icol = GLV(line,1) ; skytp[icol] = GLV(line,2) ; }} ; fclose(infu) ncc=0 for iss = -200L, 200-1 do begin sum=0. for icol=0, ecol-1 do begin ii = icol + iss if ((ii ge 0)&&(ii lt ecol)) then sum = sum + (skytp[icol] * SPX[nso].spsky[ii]) ;if ( iss eq 0) and (icol eq 1024) then print, 'nso ', nso, ' iss ', iss, ' icol ', icol, ' ncc ', ncc, ' sum ', sum endfor xcc[ncc] = iss ycc[ncc] = sum ++ncc endfor ; iss if max(ycc) eq 0 then stop ;> DEBUG arr[narr] = nsx_find_cc_peak( ncc, xcc, ycc ) ;print, "nsx_find_cc_peak -> ", arr[narr];> DEBUG ++narr endfor ccshift = cnint(( cfind_median8(narr,arr) )) print, format='("Use skyline whole pixel shift of = ",I," .")',ccshift printf,logfu, format='("Use skyline whole pixel shift of = ",I," .")',ccshift ; Skyline shift. nss=0 for nso=3,7 do begin if (nso eq 7) then ecol=IMGSKY.nc/2 else ecol=IMGSKY.nc ; Copy sky spectrum on pixel scale and apply cross correlation peak shift. */ nn=0 for icol=0L,ecol-1 do begin ii = icol + ccshift if ((ii ge 0)&&(ii lt ecol)) then begin yy[nn] = SPX[nso].spsky[ii] xx[nn] = double(icol) ++nn endif endfor ; Centroid sky lines for Sky Shift. radius=3. niter=4 fwhm_radius=4. for ii=0L,skynum[nso]-1 do begin if (skywgt[nso,ii] gt 0.) then begin xv = skywave[nso,ii] - WSC[nso].xoffinv rcol = cpolyval(WSC[nso].orderinv+1,WSC[nso].coefinv,xv) cent = nsx_centroid2( nn, xx, yy, rcol, radius, niter, 1 ) fwhm = nsx_fwhm_1D( nn, xx, yy, cent, fwhm_radius, peak, 1 ) if (cent gt -1.e+20) then begin ssnso[nss] = nso sscol[nss] = rcol sscent[nss]= cent ssresi[nss]= cent-rcol sswave[nss]= skywave[nso,ii] ssfwhm[nss]= fwhm sspeak[nss]= peak ++nss endif endif endfor endfor ; FWHM median from all orders. narr=0 for ii=0,nss-1 do begin arr[narr]=ssfwhm[ii] ++narr endfor fwhm_median = cfind_median8(narr,arr) ; Medians for different orders. */ for nso=3,7 do begin narr=0 for ii=0L,nss-1 do begin if (ssnso[ii] eq nso) then begin arr[narr]=ssresi[ii] ++narr endif endfor median = cfind_median8(narr,arr) print,format='("Echelle order ",I,": sky line shift in pixels (median)=",F9.5)',nso,median printf,logfu,format='("Echelle order ",I,": sky line shift in pixels (median)=",F9.5)',nso,median endfor ; Median from all orders. */ narr=0 for ii=0L,nss-1 do begin arr[narr]=ssresi[ii] ++narr endfor median = cfind_median8(narr,arr) ; Write skyline residuals to file. */ wrd = nsxout + IMGSKY.root + "-skyres.tbl" outfu = fopen_write(wrd) printf,outfu,"|ord| column |centroid|resid_pix|resid_ang|wavelength|shift_res|fwhm_pix| peak |fwhm_ang|" rms=0. mrms=0. ratio=0. mratio=0. for ii=0L,nss-1 do begin nso = ssnso[ii] icol= cnint(( sscent[ii] )) disp= SPX[nso].spdsp[icol] printf,outfu,format='(" ",I3," ",F8.3," ",F8.3," ",F9.4," ",F9.4," ",F10.3," ",F9.4," ",F8.3," ",F8.1," ",F8.3,"")',$ ssnso[ii], sscol[ii], sscent[ii], ssresi[ii], ssresi[ii]*disp, sswave[ii], ssresi[ii]-median,$ ssfwhm[ii], sspeak[ii], ssfwhm[ii]*ABS((disp)) rms = rms + ((ssresi[ii]-median)*(ssresi[ii]-median)) endfor fclose,outfu if (nss gt 1) then begin rms = sqrt(( rms / double(nss) )) if (rms gt 0.) then ratio = ABS((median / rms)) mrms= rms / sqrt(double(nss)) if (mrms gt 0.) then mratio = ABS((median / mrms)) skyshift = double(ccshift) + median fullratio= ABS(( skyshift / mrms )) format = '(" Skyline shift (all orders): ",F9.5,"(pixels) rms=",F9.5," [",F6.2,"] mrms=",F9.5," [",F6.2,"] {",F6.1,"} num=",I," fwhm=",F8.3)' print, format=format,skyshift,rms,ratio,mrms,mratio,fullratio,nss,fwhm_median printf, logfu, format=format,skyshift,rms,ratio,mrms,mratio,fullratio,nss,fwhm_median endif ;/* Use no skyline shift if constructing atm.abs. files. */ ;if (CONSTRUCT_ATMABS) { ; printf("WARNING: use skyline shift median of zero. \n") ; median = 0. ;} ; Apply sky shift if large enough. Apply pixel shift times dispersion. if ((nss gt 1)&&(ratio gt 0.)&&(rms gt 0.)) then begin if (fullratio gt 5.) then begin skyshift = double(ccshift) + median print,format='("Apply sky line shift of ",F9.5," pixels to all echelle orders.")',skyshift printf,logfu,format='("Apply sky line shift of ",F9.5," pixels to all echelle orders.")',skyshift for nso=3,7 do begin if (nso eq 7) then ecol=IMGSKY.nc/2 else ecol=IMGSKY.nc for icol=0L,ecol-1 do begin SPX[nso].spwav[icol] = SPX[nso].spwav[icol] - (SPX[nso].spdsp[icol] * skyshift) endfor endfor endif else begin print, "No sky line shift applied, not significant enough." printf,logfu,"No sky line shift applied, not significant enough." endelse endif ; stop ; print, "returning from nsx_compute_skyline_shift to nsx_extract_spectrum" ;> DEBUG ; return end