; docformat = 'rst' ;+ ; ; Class to manage a spectrum. ; ; Updated versions of the NASA Ames PAH IR Spectroscopic Database and ; more information can be found at: `www.astrochemistry.org/pahdb <https://www.astrochemistry.org/pahdb>`. ; ; :Examples: ; Create and destroy an ; AmesPAHdbIDLSuite_Spectrum-instance:: ; ; IDL> spectrum = OBJ_NEW('AmesPAHdbIDLSuite_Spectrum') ; IDL> spectrum->Set,data ; IDL> spectrum->Plot ; IDL> OBJ_DESTROY,spectrum ; ; :Author: ; Dr. Christiaan Boersma ; ; :Copyright: ; BSD licensed ; ; :History: ; Changes:: ; ; 04-18-2025 ; Normalize b-vector in FIT. Christiaan Boersma. ; 08-28-2024 ; Do unit conversion already in MCFIT. Christiaan Boersma. ; 08-25-2024 ; Report when maximum number of iterations reached in FIT. Christiaan ; Boersma. ; 04-09-2024 ; Normalize matrix in FIT. Christiaan Boersma. ; 11-22-2023 ; Allow specifying the tolerance and maximum number of iterations to use for ; NNLS in FIT and MCFIT and pass to result. Christiaan Boersma. ; 11-06-2023 ; Add callback to FIT and NNLS for progress monitoring. Christiaan Boersma. ; 10-19-2023 ; Don't destroy created observation in MCFIT. Christiaan Boersma. ; 04-15-2023 ; Store unperturbed observation in MCFIT. Christiaan Boersma. ; 07-08-2022 ; Pass EXTERNAL_NNLS to FIT in MCFIT. Christiaan Boersma. ; 05-18-2022 ; Use HISTOGRAM speed-up in PLOT. Christiaan Boersma. ; 04-28-2022 ; Clean up OBSERVATION in MCFIT when internally generated. ; Christiaan Boersma. ; 04-27-2022 ; Added MCFIT and corrected small typos in description of FIT. ; Christiaan Boersma. ; 08-17-2021 ; Don't try and access pointer when it is not set in DESCRIPTION. ; Christiaan Boersma. ; 05-03-2021 ; Avoid potential issues when self.uids and self.data.uid don't have ; the same ordering. Christiaan Boersma. ; 05-02-2021 ; Added NOTICE-keyword to FIT. Christiaan Boersma. ; 04-30-2021 ; Refactored FIT to make use of matrix operation. ; Christiaan Boersma. ; 09-19-2017 ; Fixed !NULL value in return for FWHM in FIT. ; 06-27-2016 ; Added EXTERNAL_NNLS-keyword to FIT. Christiaan Boersma. ; 03-23-2016 ; Removed unnecessary declaration of 'weights' in FIT. Christiaan ; Boersma. ; 03-18-2016 ; Replaced [ ] with !NULL in NNLS to avoid compiler error for IDL ; version older than 8. Christiaan Boersma. ; 11-05-2015 ; Added REGRID. Christiaan Boersma. ; 07-12-2015 ; Added method to FIT. Christiaan Boersma. ; 02-01-2015 ; First version of the file. Christiaan Boersma. ;- ;+ ; Output spectrum description. ; ; :Params: ; Str: out, optional, type="string array" ; Ouput to Str ; ; :Categories: ; INFORMATIVE ;- PRO AmesPAHdbIDLSuite_Spectrum::Description,Str COMPILE_OPT IDL2 ON_ERROR,2 self->AmesPAHdbIDLSuite_Data::Description,Str Str = [Str, STRING(FORMAT='(A-12,":",X,g-8.4,X,A-0)', "shift", self.shift, "cm!U-1!N")] Str = [Str, STRING(FORMAT='(A-12,":",X,A-0)', "profile", self.profile)] IF PTR_VALID(self.fwhm) THEN BEGIN IF SIZE(*self.fwhm, /TYPE) EQ 8 THEN BEGIN IF N_ELEMENTS((*self.fwhm)) GT 4 THEN fwhm = (*self.fwhm)[0:3].fwhm $ ELSE fwhm = (*self.fwhm).fwhm fwhm = STRJOIN(STRTRIM(STRING(FORMAT='(g-7.3)', fwhm), 2), ",") fwhm += ",..." Str = [Str, STRING(FORMAT='(A-12,":",X,A-0,X,A-0)', "FWHM", fwhm, "cm!U-1!N")] Str = [Str, STRING(FORMAT='(A-12,":",X,A-0)', "|_sectioned", "yes")] ENDIF ELSE Str = [Str, STRING(FORMAT='(A-12,":",X,g-8.4,X,A-0)', "FWHM", *self.fwhm, "cm!U-1!N")] ENDIF Str = STRJOIN(Str, "!C") IF N_PARAMS() GT 0 THEN RETURN PRINT,STRJOIN(STRSPLIT(Str, "!C", /EXTRACT, /REGEX), STRING(10B)) END ;+ ; Plot the spectrum. ; ; :Keywords: ; Wavelength: in, optional, type=int ; Whether to set the abscissa units to wavelength ; Stick: in, optional, type=int ; Whether to plot the spectrum as sticks ; Fill: in, optional, type=int ; Whether to solid-fill the spectrum ; Oplot: in, optional, type=int ; Whether to draw over a previous plot ; Legend: in, optional, type=int ; Whether to show a legend ; Color: in, optional, type=int ; Color to plot the spectrum with ; _EXTRA: in, optional, type=struct ; Required for IDL's keyword-inheritance mechanism ; ; :Categories: ; PLOTTING ;- PRO AmesPAHdbIDLSuite_Spectrum::Plot,Wavelength=Wavelength,Stick=Stick,Fill=Fill,Oplot=Oplot,Legend=Legend,Color=Color,_EXTRA=EXTRA COMPILE_OPT IDL2 ON_ERROR,2 self->AmesPAHdbIDLSuite_Plot::Setup,Oplot=Oplot,XSIZE=600,YSIZE=400 x = *self.grid nx = N_ELEMENTS(x) xunits = self.units.abscissa.str xrange = [MAX(x, MIN=xmin), xmin] IF KEYWORD_SET(Wavelength) THEN BEGIN x = 1D4 / x xrange = [MIN(x, MAX=xmax), xmax] xunits = 'wavelength [!Mm!Xm]' ENDIF IF NOT KEYWORD_SET(Oplot) THEN self->AmesPAHdbIDLSuite_Plot::Plot,REFORM(REBIN(x,nx,self.nuids),self.nuids*nx),(*self.data).intensity,Color=Color,XRANGE=xrange,XTITLE=xunits,YTITLE=self.units.ordinate.str,/NoData,Stick=Stick,_EXTRA=EXTRA IF NOT KEYWORD_SET(Color) THEN Color = 2 h = HISTOGRAM((*self.data).uid, MIN=0, REVERSE_INDICES=ri) FOR i = 0, self.nuids - 1 DO BEGIN select = ri[ri[(*self.uids)[i]]:ri[(*self.uids)[i]+1]-1] self->AmesPAHdbIDLSuite_Plot::Oplot,x,(*self.data)[select].intensity,Stick=Stick,Fill=Fill,COLOR=Color+i ENDFOR IF SIZE(Legend, /TYPE) EQ 0 THEN Legend = 1 IF Legend THEN BEGIN self->Description,outs self->AmesPAHdbIDLSuite_Plot::Legend,outs ENDIF self->AmesPAHdbIDLSuite_Plot::Restore END ;+ ; Write the spectrum to file as an IPAC-table. ; ; :Params: ; Filename: in, optional, type=string ; Output filename ; ; :Categories: ; OUTPUT ;- PRO AmesPAHdbIDLSuite_Spectrum::Write,Filename COMPILE_OPT IDL2 ON_ERROR,2 IF N_PARAMS() LT 1 THEN Filename = OBJ_CLASS(self) + '.tbl' timestamp = SYSTIME() hdr = [] FXADDPAR,hdr,"DATE",timestamp," Date this file was generated" FXADDPAR,hdr,"ORIGIN","NASA Ames Research Center"," Organization generating this file" FXADDPAR,hdr,"CREATOR",STRING(FORMAT='("IDL",X,A0,X,"on",X,A0)', !VERSION.RELEASE, !VERSION.OS_NAME)," Software used to create this file" FXADDPAR,hdr,"SOFTWARE","AmesPAHdbIDLSuite"," Program used to create this file" FXADDPAR,hdr,"AUTHOR","Dr. C. Boersma"," Author of the program" FXADDPAR,hdr,"TYPE",OBJ_CLASS(self)," AmesPAHdbIDLSuite data type" self->Description,description comments = STRSPLIT(description, "!C", /EXTRACT, /REGEX, COUNT=ncomments) FOR i = 0L, ncomments - 1 DO FXADDPAR,hdr,"COMMENT",comments[i] IF self.units.abscissa.str THEN $ abscissa = STREGEX(self.units.abscissa.str, '(.*) \[(.*)\]', /SUBEXPR, /EXTRACT) $ ELSE $ abscissa = ['', 'abscissa', ''] IF self.units.ordinate.str THEN $ ordinate = STREGEX(self.units.ordinate.str, '(.*) \[(.*)\]', /SUBEXPR, /EXTRACT) $ ELSE $ ordinate = ['', 'ordinate', ''] half_abscissa_len = STRLEN(abscissa[1]) / 2 half_ordinate_len = STRLEN(ordinate[1]) / 2 fmt1 = '("|",A' + STRING(FORMAT='(I0)', 12 + half_abscissa_len) + ',' + STRING(FORMAT='(I0)', 13 - half_abscissa_len) + 'X,' + $ '"|",A' + STRING(FORMAT='(I0)', 12 + half_ordinate_len) + ',' + STRING(FORMAT='(I0)', 13 - half_ordinate_len) + 'X,' + $ '"|",A' + STRING(FORMAT='(I0)', 6 + 3) + ',' + STRING(FORMAT='(I0)', 6 - 3) + 'X,' + $ '"|")' fmt2 = '("|",A' + STRING(FORMAT='(I0)', 12 + 3) + ',' + STRING(FORMAT='(I0)', 13 - 3) + 'X,' + $ '"|",A' + STRING(FORMAT='(I0)', 12 + 3) + ',' + STRING(FORMAT='(I0)', 13 - 3) + 'X,' + $ '"|",A' + STRING(FORMAT='(I0)', 6 + 3) + ',' + STRING(FORMAT='(I0)', 6 - 3) + 'X,' + $ '"|")' half_abscissa_len = STRLEN(abscissa[2]) / 2 half_ordinate_len = STRLEN(ordinate[2]) / 2 fmt3 = '("|",A' + STRING(FORMAT='(I0)', 12 + half_abscissa_len) + ',' + STRING(FORMAT='(I0)', 13 - half_abscissa_len) + 'X,' + $ '"|",A' + STRING(FORMAT='(I0)', 12 + half_ordinate_len) + ',' + STRING(FORMAT='(I0)', 13 - half_ordinate_len) + 'X,' + $ '"|",A' + STRING(FORMAT='(I0)', 6 + 3) + ',' + STRING(FORMAT='(I0)', 6 - 3) + 'X,' + $ '"|")' cols = [STRING(FORMAT=fmt1,STRUPCASE(abscissa[1]),STRUPCASE(ordinate[1]),'UID'), $ STRING(FORMAT=fmt2,"double","double","int"), $ STRING(FORMAT=fmt3,abscissa[2],ordinate[2],"")] n = N_ELEMENTS(*self.grid) intensities = REFORM((*self.data).intensity, n, self.nuids) srt = SORT(*self.uids) intensities = intensities[*, srt] OPENW,funit,Filename,/GET_LUN PRINTF,funit,FORMAT='("\",A0)',hdr[0:WHERE(STRPOS(hdr, 'END') EQ 0)] PRINTF,funit,STRJOIN(cols, STRING( 10B )) n = N_ELEMENTS(*self.grid) FOR i = 0L, self.nuids - 1L DO BEGIN FOR j = 0L, n - 1L DO PRINTF,funit,FORMAT='(X,F25.6,X,F25.6,X,I)',(*self.grid)[j],intensities[j,i],(*self.uids)[i] ENDFOR CLOSE,funit FREE_LUN,funit PRINT PRINT,"=========================================================" PRINT," WRITTEN IPAC TABLE: ", Filename PRINT,"=========================================================" PRINT END ;+ ; Implementation of the non-negative least-squares algorithm. ; ; :Params: ; A: in, required, type="double array (2D)" ; Matrix ; b: in, required, type="double array (1D)" ; Vector ; tol: in, required, type=double ; Tolerance ; max_iter: in, required, type=long ; Maximum number of iterations ; ; :Keywords: ; callback: in, optional, type=string ; Procedure called on every iteration ; ; :Categories: ; FITTING ; ; :Private: ;- PRO AmesPAHdbIDLSuite_Spectrum::NNLS,A,b,tol,max_iter,callback=callback COMPILE_OPT IDL2 ON_ERROR, 2 A = DOUBLE(A) b = DOUBLE(b) tol = DOUBLE(tol) max_iter = LONG(max_iter) s_struct = SIZE(A, /STRUCTURE) m = s_struct.dimensions[0] ; n = s_struct.dimensions[1] w = MAKE_ARRAY(m, VALUE=-1D) x = DBLARR(m) P = !NULL Z = INDGEN(m) k = 0L WHILE (Z NE !NULL) AND (k LT max_iter) DO BEGIN IF WHERE(-w[Z] GT tol, /NULL) EQ !NULL THEN BREAK k += 1L tapbp = A#b IF P EQ !NULL THEN BEGIN w = -tapbp ENDIF ELSE BEGIN tapap = A#TRANSPOSE(A[P, *]) w = -tapbp + tapap#x[P] ENDELSE IF WHERE(-w[Z] GT tol, /NULL) EQ !NULL THEN BREAK wi = MAX(-w[Z], i) P = [P, Z[i]] nz = N_ELEMENTS(Z) IF nz EQ 1 THEN BEGIN Z = !NULL ENDIF ELSE BEGIN Z = Z[WHERE(HISTOGRAM([i], MIN=0, MAX=nz-1) EQ 0)] ENDELSE WHILE 1 DO BEGIN zz = DBLARR(m) zz[P] = LA_LEAST_SQUARES(A[P,*],b) IF WHERE(zz[P] LE 0, /NULL) EQ !NULL THEN BEGIN x = zz BREAK ENDIF P_ZN = P[WHERE(zz[P] LE tol, /NULL)] alpha = MIN(x[P_ZN]/(x[P_ZN] - zz[P_ZN])) x[P] = x[P] + alpha*(zz[P]-x[P]) temp = WHERE(ABS(x[P]) LE tol, /NULL) Z = [Z, P[temp]] np = N_ELEMENTS(P) IF np EQ N_ELEMENTS(temp) THEN BEGIN P = !NULL ENDIF ELSE IF temp NE !NULL THEN BEGIN P = P[WHERE(HISTOGRAM(temp, MIN=0, MAX=np-1) EQ 0)] ENDIF ENDWHILE IF KEYWORD_SET(callback) THEN CALL_PROCEDURE,callback,A,x ENDWHILE b=x max_iter=k END ;+ ; Perform a spectroscopic fit. ; ; :Returns: ; AmesPAHdbIDLSuite_Fitted_Spectrum ; ; :Params: ; observation: in, required, type="double array (1D) or AmesPAHdbIDLSuite_Observation" ; Observed spectrum ; error: in, optional, type="double array (1D)" ; Uncertainties associated with observation ; ; :Keywords: ; TOLERANCE_NNLS: in, optional, type=double ; Tollerance ; MAXITER_NNLS: in, optional, type=long ; Maximum number of iterations ; CALLBACK_NNLS: in, optional, type=boolean ; Callback ; EXTERNAL_NNLS: in, optional, type=int ; Whether to use an externally defined NNLS-routine ; Notice: in, optional, type=int, default=1 ; Whether to show notices ; ; :Categories: ; FITTING ;- FUNCTION AmesPAHdbIDLSuite_Spectrum::Fit,observation,error,TOLERANCE_NNLS=tolerance_nnls,MAXITER_NNLS=maxiter_nnls,CALLBACK_NNLS=callback_nnls,EXTERNAL_NNLS=external_nnls,Notice=Notice COMPILE_OPT IDL2 ON_ERROR,2 IF SIZE(Notice, /TYPE) EQ 0 THEN Notice = 1 IF SIZE(TOLERANCE_NNLS, /TYPE) EQ 0 THEN tolerance_nnls = 1D-16 ftol = tolerance_nnls IF SIZE(MAXITER_NNLS, /TYPE) EQ 0 THEN maxiter_nnls = 128L type = SIZE(observation, /STRUCTURE) has_error = 0 IF type.type_name EQ 'OBJREF' THEN BEGIN IF OBJ_CLASS(observation) EQ 'AMESPAHDBIDLSUITE_OBSERVATION' THEN BEGIN observation->AbscissaUnitsTo,1,Notice=Notice observation_s = observation->get() IF NOT ARRAY_EQUAL(*self.grid, observation_s.data.x) THEN BEGIN PRINT PRINT,"=========================================================" PRINT," DATA AND OBSERVATION GRIDS ARE NOT THE SAME " PRINT,"=========================================================" PRINT RETURN,OBJ_NEW() ENDIF IF TOTAL(observation_s.data.ystdev) GT 0 THEN has_error = 1 ENDIF ELSE BEGIN PRINT PRINT,"=========================================================" PRINT," OBJECT SHOULD BE AN AMESPAHDBIDLSUITE_OBSERVATION: "+type.type_name PRINT,"=========================================================" PRINT self.state = 0 RETURN,OBJ_NEW() ENDELSE ENDIF ELSE BEGIN IF N_PARAMS() GT 1 THEN BEGIN tmp = OBJ_NEW('AmesPAHdbIDLSuite_Observation', $ X=*self.grid, $ Y=observation, $ ErrY=error) has_error = 1 ENDIF ELSE tmp = OBJ_NEW('AmesPAHdbIDLSuite_Observation', $ X=*self.grid, $ Y=observation) observation_s = tmp->get() OBJ_DESTROY,tmp ENDELSE ny = N_ELEMENTS(observation_s.data.y) matrix = TRANSPOSE(REFORM((*self.data).intensity, ny, self.nuids)) m = matrix b = observation_s.data.y - observation_s.data.continuum IF has_error THEN BEGIN IF Notice THEN BEGIN PRINT PRINT,"=========================================================" PRINT," DOING NNLC " PRINT,"=========================================================" PRINT ENDIF b /= observation_s.data.ystdev m /= TRANSPOSE(REBIN(observation_s.data.ystdev, ny, self.nuids)) method = 'NNLC' ENDIF ELSE BEGIN IF Notice THEN BEGIN PRINT PRINT,"=========================================================" PRINT," DOING NNLS " PRINT,"=========================================================" PRINT ENDIF method = 'NNLS' ENDELSE READS,!VERSION.RELEASE,idl_version b_scl = MAX(b) b /= b_scl m_scl = MAX(m) m /= m_scl IF idl_version GE 8.0 AND NOT KEYWORD_SET(EXTERNAL_NNLS) THEN BEGIN maxiter = maxiter_nnls self->NNLS,m,b,ftol,maxiter_nnls,callback=callback_nnls IF maxiter EQ maxiter_nnls AND Notice THEN BEGIN PRINT PRINT,"=========================================================" PRINT," MAXIMUM NUMBER OF ITERATIONS REACHED: "+STRTRIM(maxiter,2) PRINT,"=========================================================" PRINT ENDIF weights = b ENDIF ELSE IF NOT KEYWORD_SET(EXTERNAL_NNLS) THEN BEGIN PRINT PRINT,"=========================================================" PRINT," FIT REQUIRES IDL VERSION 8.0 OR HIGHER: "+!VERSION.RELEASE PRINT,"=========================================================" PRINT self.state = 0 RETURN,OBJ_NEW() ENDIF ELSE BEGIN IF Notice THEN BEGIN PRINT PRINT,"=========================================================" PRINT," USING EXTERNAL NNLS " PRINT,"=========================================================" PRINT ENDIF weights = DBLARR(self.nuids, /NOZERO) enorm = 0D w = DBLARR(self.nuids) indx = LONARR(self.nuids) mode = 0 NNLS,m,ny,self.nuids,b,weights,enorm,w,indx,mode ENDELSE valid = WHERE(weights GT 0, nvalid) IF nvalid EQ 0 THEN BEGIN PRINT PRINT,"=========================================================" PRINT," UNABLE TO FIND SOLUTION " PRINT,"=========================================================" PRINT self.state = 0 RETURN,OBJ_NEW() ENDIF uids = (*self.data)[UNIQ((*self.data).uid)].uid _weights = REPLICATE({AmesPAHdbIDLSuite_Weights_S, $ uid:0L, $ weight:0D}, nvalid) _weights.uid = uids[valid] _weights.weight = weights[valid] * b_scl / m_scl data = REPLICATE({AmesPAHdbIDLSuite_Fitted_S, $ intensity:0D, $ uid:0L}, ny * nvalid) FOR i = 0, nvalid - 1 DO BEGIN data[i*ny:(i+1)*ny-1].uid = _weights[i].uid data[i*ny:(i+1)*ny-1].intensity = _weights[i].weight * REFORM(matrix[valid[i], *], /OVERWRITE) ENDFOR IF Notice THEN BEGIN PRINT PRINT,"=========================================================" PRINT," NOTICE: PLEASE TAKE CONSIDERABLE CARE WHEN INTERPRETING " PRINT," THESE RESULTS AND PUTTING THEM IN AN ASTRONOMICAL " PRINT," CONTEXT. THERE ARE MANY SUBTLETIES THAT NEED TO BE TAKEN" PRINT," INTO ACCOUNT, RANGING FROM PAH SIZE, INCLUSION OF " PRINT," HETEROATOMS, ETC. TO DETAILS OF THE APPLIED EMISSION " PRINT," MODEL, BEFORE ANY THOROUGH ASSESSMENT CAN BE MADE. " PRINT,"=========================================================" PRINT ENDIF RETURN,OBJ_NEW('AmesPAHdbIDLSuite_Fitted_Spectrum', $ Type=self.type, $ Version=self.version, $ Data=data, $ PAHdb=self.database, $ Uids=_weights.uid, $ Model=*self.model, $ Units=self.units, $ Shift=self.shift, $ Grid=*self.grid, $ Profile=self.profile, $ FWHM=self.fwhm NE !NULL ? *self.fwhm : !NULL, $ Observation=observation_s, $ Weights=_weights, $ Method=method, $ Tolerance=ftol, $ Iterations=maxiter_nnls) END ;+ ; Perform spectroscopic fits using a Monte Carlo approach. ; ; :Returns: ; AmesPAHdbIDLSuite_MCFitted_Spectrum ; ; :Params: ; observation: in, required, type="double array (1D) or AmesPAHdbIDLSuite_Observation" ; Observed spectrum ; error: in, required, type="double array (1D)" ; Uncertainties associated with observation ; samples: in, required, type=int ; Number of Monte Carlo samples ; ; :Keywords: ; Uniform: in, optional, type=int ; Whether to use a uniform rather than a normal distribution to permutate the errors. ; TOLERANCE_NNLS: in, optional, type=double ; Tollerance ; MAXITER_NNLS: in, optional, type=long ; Maximum number of iterations ; EXTERNAL_NNLS: in, optional, type=int ; Whether to use an externally defined NNLS-routine. ; ; :Categories: ; FITTING ;- FUNCTION AmesPAHdbIDLSuite_Spectrum::MCFit,observation,error,samples,TOLERANCE_NNLS=tolerance_nnls,MAXITER_NNLS=maxiter_nnls,EXTERNAL_NNLS=external_nnls,Uniform=Uniform COMPILE_OPT IDL2 ON_ERROR,2 type = SIZE(observation, /STRUCTURE) IF type.type_name EQ 'OBJREF' THEN BEGIN samples = error obs = observation obs->AbscissaUnitsTo,1 obs_s = obs->Get() y = obs_s.data.y ystdev = obs_s.data.ystdev ENDIF ELSE BEGIN IF N_PARAMS() GT 1 THEN $ obs = OBJ_NEW('AmesPAHdbIDLSuite_Observation', $ X=*self.grid, $ Y=observation, $ ErrY=error) $ ELSE $ obs = OBJ_NEW('AmesPAHdbIDLSuite_Observation', $ X=*self.grid, $ Y=observation) y = observation ystdev = error ENDELSE ny = N_ELEMENTS(y) obj = OBJARR(samples) PRINT PRINT,"=========================================================" PRINT," DOING MONTE CARLO " PRINT,"=========================================================" PRINT PRINT PRINT,"=========================================================" IF KEYWORD_SET(Uniform) THEN $ PRINT," DRAWING FROM UNIFORM DISTRIBUTION " $ ELSE $ PRINT," DRAWING FROM NORMAL DISTRIBUTION " PRINT,"=========================================================" PRINT PRINT PRINT,"=========================================================" FOR i = 0L, samples - 1L DO BEGIN PRINT,FORMAT='("' + STRING(13B) + 'mc:",X,I5,"/",I5,$)',i+1L,samples obs->Set,Y=y + ystdev * (NOT KEYWORD_SET(Uniform) ? RANDOMU(seed, ny, /DOUBLE, /NORMAL) $ : (2D * RANDOMU(seed, ny, /DOUBLE, /UNIFORM) - 1D)) obj[i] = self->Fit(obs, TOLERANCE_NNLS=tolerance_nnls, MAXITER_NNLS=maxiter_nnls, EXTERNAL_NNLS=external_nnls, Notice=0B) ENDFOR PRINT PRINT,"=========================================================" obs->Set,Y=y RETURN,OBJ_NEW('AmesPAHdbIDLSuite_MCFitted_Spectrum', $ Type=self.type, $ Obj=obj, $ Distribution=KEYWORD_SET(Uniform) ? 'uniform' : 'normal', $ Observation=obs) END ;+ ; Resample the spectrum onto a provided grid. ; ; :Params: ; grid: in, required, type="double array (1D)" ; The grid ; ; :Categories: ; MANIPULATE ;- PRO AmesPAHdbIDLSuite_Spectrum::Resample,grid COMPILE_OPT IDL2 ON_ERROR,2 IF NOT PTR_VALID(self.data) THEN BEGIN PRINT PRINT,"=========================================================" PRINT," NO DATA " PRINT,"=========================================================" PRINT self.state = 0 RETURN ENDIF ngrid1 = N_ELEMENTS(grid) Data = REPLICATE({AmesPAHdbSpectrum_S, $ intensity:0D, $ uid:0L}, ngrid1 * self.nuids) ngrid2 = N_ELEMENTS(*self.grid) uids = (*self.data)[UNIQ((*self.data).uid)].uid i = 0 FOR j = 0, ngrid1 - 1 DO BEGIN idx = i WHILE i + 1 LT ngrid2 - 1 DO BEGIN IF (*self.grid)[i+1] GT grid[j] THEN BREAK idx = [idx, ++i] ENDWHILE FOR k = 0, self.nuids - 1 DO BEGIN Data[k*ngrid1+j].uid = uids[k] IF N_ELEMENTS(idx) EQ 1 THEN BEGIN IF (*self.grid)[i] GE grid[0] AND (*self.grid)[i] LE grid[ngrid1-1] THEN Data[k*ngrid1+j].intensity = (*self.data)[k*ngrid2+i].intensity ENDIF ELSE Data[k*ngrid1+j].intensity = MEAN((*self.data)[k*ngrid2+idx].intensity) ENDFOR ENDFOR PTR_FREE,self.grid self.grid = PTR_NEW(grid) PTR_FREE,self.data self.data = PTR_NEW(Data) END ;+ ; Co-adds the spectra ; ; :Returns: ; AmesPAHdb_Coadded_Spectrum ; ; :Keywords: ; Weights: in, optional, type=struct ; Use the provided weights when co-adding ; Average: in, optional, type=int ; Take the average ; ; :Categories: ; CALCULATE ;- FUNCTION AmesPAHdbIDLSuite_Spectrum::Coadd,Weights=weights,Average=Average COMPILE_OPT IDL2 ON_ERROR,2 select1 = WHERE((*self.data).uid EQ (*self.uids)[0], nselect1) data = REPLICATE({AmesPAHdbIDLSuite_Coadd_S, intensity:0D}, nselect1) data.intensity = (*self.data)[select1].intensity IF KEYWORD_SET(Weights) THEN BEGIN select2 = WHERE(Weights.uid EQ (*self.uids)[0], nselect2) IF nselect2 EQ 0 THEN BEGIN PRINT PRINT,"=========================================================" PRINT," NO WEIGHT DEFINED FOR UID: "+STRING(FORMAT='(I-0)',(*self.uids)[0]) PRINT,"=========================================================" PRINT self.state = 0 RETURN,0 ENDIF data.intensity *= Weights[select2].weight ENDIF FOR i = 1, self.nuids - 1 DO BEGIN select1 = WHERE((*self.data).uid EQ (*self.uids)[i]) IF KEYWORD_SET(Weights) THEN BEGIN select2 = WHERE(Weights.uid EQ (*self.uids)[i], nselect2) IF nselect2 EQ 0 THEN BEGIN PRINT PRINT,"=========================================================" PRINT," NO WEIGHT DEFINED FOR UID: "+STRING(FORMAT='(I-0)',(*self.uids)[i]) PRINT,"=========================================================" PRINT self.state = 0 RETURN,0 ENDIF data.intensity += Weights[select2].weight * (*self.data)[select1].intensity ENDIF ELSE data.intensity += (*self.data)[select1].intensity ENDFOR IF KEYWORD_SET(Average) THEN data.intensity /= self.nuids RETURN,OBJ_NEW('AmesPAHdbIDLSuite_Coadded_Spectrum', $ Type=self.type, $ Version=self.version, $ Data=data, $ PAHdb=self.database, $ Uids=*self.uids, $ Model=*self.model, $ Units=self.units, $ Shift=self.shift, $ Grid=*self.grid, $ Profile=self.profile, $ FWHM=*self.fwhm, $ Weights=Weights, $ Averaged=KEYWORD_SET(Average)) END ;+ ; Retrieves the AmesPAHdbIDLSuite_Spectrum representation in a ; structure. ; ; :Returns: ; Structure ; ; :Categories: ; SET/GET ;- FUNCTION AmesPAHdbIDLSuite_Spectrum::Get COMPILE_OPT IDL2 ON_ERROR,2 IF NOT PTR_VALID(self.data) THEN RETURN, 0 struct = self->AmesPAHdbIDLSuite_Data::Get() struct.type = OBJ_CLASS(self)+'_S' RETURN,CREATE_STRUCT(struct, 'shift', self.shift, 'grid', *self.grid, 'profile', self.profile, 'fwhm', *self.fwhm) END ;+ ; Populates the AmesPAHdbIDLSuite_Spectrum-instance. ; ; :Params: ; Struct: in, optional, type=struct ; Data structure ; ; :Keywords: ; Type: in, optional, type=string ; Type of Data ; Version: in, optional, type=string ; Versioning information ; Data: in, optional, type=struct ; Data structure ; PAHdb: in, optional, type=pointer ; Pointer to parsed database file ; Uids: in, optional, type="long array (1D)" ; UIDs in Data ; Model: in, optional, type=string ; References ; Units: in, optional, type="AmesPAHdb_Units_S struct" ; Units ; Shift: in, optional, type=float ; Shift ; Grid: in, optional, type="float array" ; Grid ; Profile: in, optional, type=string ; Profile ; FWHM: in, optional, type=float ; FWHM ; ; :Categories: ; SET/GET ;- PRO AmesPAHdbIDLSuite_Spectrum::Set,Struct,Type=Type,Version=Version,Data=Data,PAHdb=PAHdb,Uids=Uids,Model=Model,Units=Units,Shift=Shift,Grid=Grid,Profile=Profile,FWHM=FWHM COMPILE_OPT IDL2 ON_ERROR,2 IF N_PARAMS() GT 0 THEN BEGIN tag = WHERE(TAG_NAMES(Struct) EQ 'TYPE', ntype) IF ntype EQ 1 THEN BEGIN IF Struct.(tag) EQ OBJ_CLASS(self)+'_S' THEN BEGIN IF NOT KEYWORD_SET(Shift) THEN self.shift = Struct.shift IF NOT KEYWORD_SET(Grid) THEN BEGIN IF PTR_VALID(self.grid) THEN PTR_FREE,self.grid self.grid = PTR_NEW(Struct.grid) ENDIF IF NOT KEYWORD_SET(Profile) THEN self.profile = Struct.profile IF NOT KEYWORD_SET(FWHM) THEN BEGIN IF PTR_VALID(self.fwhm) THEN PTR_FREE,self.fwhm self.fwhm = PTR_NEW(Struct.fwhm) ENDIF s = Struct s.type = 'AMESPAHDBIDLSUITE_Data_S' self->AmesPAHdbIDLSuite_Data::Set,Struct,Type=Type,Version=Version,Data=Data,PAHdb=PAHdb,Uids=Uids,Model=Model,Units=Units ENDIF ENDIF ENDIF ELSE self->AmesPAHdbIDLSuite_Data::Set,Type=Type,Version=Version,Data=Data,PAHdb=PAHdb,Uids=Uids,Model=Model,Units=Units IF KEYWORD_SET(Shift) THEN self.shift = Shift IF KEYWORD_SET(Grid) THEN BEGIN IF PTR_VALID(self.grid) THEN PTR_FREE,self.grid self.grid = PTR_NEW(Grid) ENDIF IF KEYWORD_SET(Profile) THEN self.profile = Profile IF KEYWORD_SET(FWHM) THEN BEGIN IF PTR_VALID(self.fwhm) THEN PTR_FREE,self.fwhm self.fwhm = PTR_NEW(FWHM) ENDIF END ;+ ; Retrieves the abscissa values. ; ; :Returns: ; double array (1D) ; ; :Categories: ; SET/GET ;- FUNCTION AmesPAHdbIDLSuite_Spectrum::GetGrid COMPILE_OPT IDL2 ON_ERROR,2 IF PTR_VALID(self.grid) THEN RETURN, *self.grid RETURN,0 END ;+ ; Clean-up an AmesPAHdbIDLSuite_Spectrum-instance ; ; :Categories: ; CLASS ; ; :Private: ;- PRO AmesPAHdbIDLSuite_Spectrum::Cleanup COMPILE_OPT IDL2 ON_ERROR,2 self->AmesPAHdbIDLSuite_Plot::Cleanup self->AmesPAHdbIDLSuite_Data::Cleanup IF PTR_VALID(self.grid) THEN PTR_FREE,self.grid IF PTR_VALID(self.fwhm) THEN PTR_FREE,self.fwhm END ;+ ; Create an AmesPAHdbIDLSuite_Spectrum-instance ; ; :Returns: ; AmesPAHdbIDLSuite_Spectrum-instance ; ; :Params: ; Struct: in, optional, type=struct ; Data structure ; ; :Keywords: ; Type: in, optional, type=string ; Type of Data ; Version: in, optional, type=string ; Versioning information ; Data: in, optional, type=struct ; Data structure ; PAHdb: in, optional, type=pointer ; Pointer to parsed database file ; Uids: in, optional, type="long array (1D)" ; UIDs in Data ; Model: in, optional, type=string ; References ; Units: in, optional, type="AmesPAHdb_Units_S struct" ; Units ; Shift: in, optional, type=float ; Shift ; Grid: in, optional, type="float array" ; Grid ; Profile: in, optional, type=string ; Profile ; FWHM: in, optional, type=float ; FWHM ; ; :Categories: ; CLASS ;- FUNCTION AmesPAHdbIDLSuite_Spectrum::Init,Struct,Type=Type,Version=Version,Data=Data,PAHdb=PAHdb,Uids=Uids,Model=Model,Units=Units,Shift=Shift,Grid=Grid,Profile=Profile,FWHM=FWHM COMPILE_OPT IDL2 ON_ERROR,2 self.state = self->AmesPAHdbIDLSuite_Plot::Init() IF self.state EQ 1 THEN BEGIN IF N_PARAMS() GT 0 THEN self->Set,Struct,Type=Type,Version=Version,Data=Data,PAHdb=PAHdb,Uids=Uids,Model=Model,Units=Units,Shift=Shift,Grid=Grid,Profile=Profile,FWHM=FWHM $ ELSE self->Set,Type=Type,Version=Version,Data=Data,PAHdb=PAHdb,Uids=Uids,Model=Model,Units=Units,Shift=Shift,Grid=Grid,Profile=Profile,FWHM=FWHM ENDIF RETURN,self.state END ;+ ; Defines the AmesPAHdbIDLSuite_Spectrum Class ; ; :Fields: ; shift: type=double ; Applied band shift ; grid: type=pointer ; Pointer to abscissa values ; profile: type=string ; Applied band profile ; fwhm: type=pointer ; Pointer to the FWHM of the applied band profile ; ; :Categories: ; CLASS ; ; :Private: ;- PRO AmesPAHdbIDLSuite_Spectrum__DEFINE COMPILE_OPT IDL2 ON_ERROR,2 void = {AmesPAHdbIDLSuite_Spectrum, $ INHERITS AmesPAHdbIDLSuite_Plot, $ INHERITS AmesPAHdbIDLSuite_Data, $ shift:0D, $ grid:PTR_NEW(), $ profile:'', $ fwhm:PTR_NEW()} END ; END OF amespahdbidlsuite_spectrum__define.pro