@pikaia.pro ; TODO: insert cos(i) in hankel transform? ; ////////////////////////////////////////////////////////////////////////////////////////// ; Read the parameter file ; ////////////////////////////////////////////////////////////////////////////////////////// function read_param_file, param_file,ndust, geometry, fit_individual_sizes, n_zones common obs_params, dist, rstar, tstar, nwav, nwav_all, owavl, oflux, ouflux, fit_range, source_name, bll, blpa, rsub, save_chain,lun,fit_incl readcol, param_file, names, values, format='(A,A)', delimiter='=' ;['Inner radius of the rim [AU]'],$ ;['Inner radius of the disk [AU]'],$ param_name = [ ['Outer radius of the disk [AU]'],$ ;['Number of points in the rim'],$ ['Number of points in the disk'],$ ;['Number of parameters to be fitted'],$ ['Number of individuals in the pop.'],$ ['Number of generations'],$ ['Number of singificant digits'],$ ['Crossover probability'],$ ['Mutation mode'],$ ['Initial mutation rate'],$ ['Minimum mutation rate'],$ ['Maximum mutation rate'],$ ['Relative fitness differential'],$ ['Reproduction plan'],$ ['Elitism flag'],$ ['Printed output (none/min/verbose)'] ] ;['Bounds for Trim'],$ bound_name = [ ['Bounds for Tdust'],$ ['Bounds for Tcont'],$ ;['Bounds for qrim'],$ ['Bounds for pdust'],$ ['Bounds for pcont'],$ ['Bounds for R1'], $ ['Bounds for R2']] ;['Bounds for r_gap_in'],$ ;['Bounds for r_gap_out']] if n_zones eq 3 then begin bound_name=[[bound_name],[['Bounds for R3']]] endif ;if fit_individual_sizes eq 1 then begin bound_name=[[bound_name],[['Bounds for q_surfdens1']]] bound_name=[[bound_name],[['Bounds for q_surfdens2']]] ;endif if n_zones eq 3 then begin bound_name=[[bound_name],[['Bounds for q_surfdens3']]] endif if fit_incl eq 1 then begin bound_name=[[bound_name],[['Bounds for cos_i']]] bound_name=[[bound_name],[['Bounds for PA']]] endif param_val = fltarr(n_elements(param_name))*0.0 - 1 if fit_individual_sizes eq 1 then begin bound_val = fltarr(2,n_elements(bound_name)-n_zones + n_zones*ndust)*0.0 - 1 endif else begin bound_val = fltarr(2,n_elements(bound_name))*0.0 - 1 endelse for j = 0, n_elements(names)-1 do begin ;print,j,names[j],n_elements(names) ; Scan the next line from the file (names[j] and values[j]) idl = strpos(names[j], '<') idr = strpos(names[j], '>') text = strmid(names[j], idl+1, (idr-idl)-1) ; Check if it was a parameter for i=0, n_elements(param_name)-1 do begin if text eq param_name[i] then begin param_val[i] = float(values[j]) print,text,i,param_val[i] endif endfor ; Check if it was a boundary for i=0, n_elements(bound_name)-1 do begin if text eq bound_name[i] then begin dum3 = strsplit(values[j], /extract) if fit_individual_sizes eq 1 then begin if text eq 'Bounds for q_surfdens1' then begin for k=0,ndust-1 do begin bound_val[0,-ndust+k] = dum3[0] bound_val[1,-ndust+k] = dum3[1] print,text,-ndust+k,' ',dum3[0],' ',dum3[1] endfor endif else begin if text eq 'Bounds for q_surfdens2' then begin for k=0,ndust-1 do begin bound_val[0,-2*ndust+k] = dum3[0] bound_val[1,-2*ndust+k] = dum3[1] print,text,-2*ndust+k,' ',dum3[0],' ',dum3[1] endfor endif else begin if text eq 'Bounds for q_surfdens3' then begin for k=0,ndust-1 do begin bound_val[0,-3*ndust+k] = dum3[0] bound_val[1,-3*ndust+k] = dum3[1] print,text,-3*ndust+k,' ',dum3[0],' ',dum3[1] endfor endif else begin bound_val[0,i] = dum3[0] bound_val[1,i] = dum3[1] endelse endelse print,text,i,' ',dum3[0],' ',dum3[1] endelse endif else begin print,text,i,' ',dum3[0],' ',dum3[1] bound_val[0,i] = dum3[0] bound_val[1,i] = dum3[1] endelse endif endfor endfor print,'b',bound_val ctrl = param_val[2:13] ;(6:17) ;rin = param_val(1) ;!!! ;rrim = param_val(0) ;!!! rout = param_val(0) ;nr_rim = param_val(3) nr_disk = param_val(1) nvar = n_elements(bound_val[0,*]) ;if fit_individual_sizes eq 1 then begin ; nvar = param_val(5) + 2.0*ndust ;endif else begin ; nvar = param_val(5) + 2 ;endelse print,'Number of free parameters: ',nvar print,'Number of dust species: ',ndust result = {ctrl:ctrl, nvar:nvar, rout:rout, nr_disk:nr_disk, pbound:bound_val} return, result end ; ////////////////////////////////////////////////////////////////////////////////////////// ; Initialize the fitting (read the opacities, PAH-templates, fitting parameter ranges, etc.) ; ////////////////////////////////////////////////////////////////////////////////////////// function read_opac, opac_fname, do_interpolate common opacities, kappa_lst, ndust, qval_fname, pah_fname, npah, pah_prof,myst_ix common obs_params, dist, rstar, tstar, nwav, nwav_all, owavl, oflux, ouflux, fit_range, source_name, bll, blpa, rsub, save_chain,lun,fit_incl print,'Reading opacity inp file: ',+opac_fname openr, 1, opac_fname qval_path = '' readf, 1, qval_path nall = 0 readf, 1, nall, format='(F)' ;print,'a',qval_path,nall qval_fname = '' for i=0, nall-1 do begin dum = '' readf, 1, dum ; Each entry consists of a number (1/0) and the filename ; The file will only be used if the number != 0 sdum = strsplit(dum, /extract) dum1 = fix(sdum[0]) dum2 = sdum[1] if dum1 eq 1 then qval_fname = [qval_fname, dum2] endfor qval_fname = qval_fname(1:n_elements(qval_fname)-1) ndust = n_elements(qval_fname) print,'b',ndust ; ; Skip any number of empty lines ; dum1 = '' dum2 = '' while dum1 eq '' do readf, 1, dum1 ; ; Read the PAH template filenames, dum1 is their directory ; readf, 1, dum2 pah_path = dum1 npah = fix(dum2) pah_fname = '' for i=0, npah-1 do begin dum = '' readf, 1, dum ; Each entry consists of a number (1/0) and a filename ; The fiel is only used if number != 0 sdum = strsplit(dum, /extract) dum1 = fix(sdum[0]) dum2 = sdum[1] if dum1 eq 1 then pah_fname = [pah_fname, dum2] endfor if n_elements(pah_fname) eq 1 and pah_fname(0) eq '' then begin npah = 0 pah_fname = '' endif else begin pah_fname = pah_fname(1:n_elements(pah_fname)-1) npah = n_elements(pah_fname) endelse close, 1 free_lun, 1 ; ; Compatibility check with the Fortran90 version ; ii = strpos(qval_path, "'") if ii(0) ge 0 then begin qval_path = strsplit(qval_path, "'", /extract) if n_elements(qval_path) ne 1 then qval_path = qval_path(1) endif ii = strpos(pah_path, "'") if ii(0) ge 0 then begin pah_path = strsplit(pah_path, "'", /extract) if n_elements(pah_path) ne 1 then pah_path = pah_path(1) endif ; ; Read the opacities ; print,'Read the opacities' kappa_lst = list() myst_ix = -1 for i = 0,n_elements(nwav)-1 do begin if do_interpolate then begin kappa = fltarr(nwav[i], ndust) endif else begin readcol, qval_path[0] + qval_fname[0], dum1, dum2, skipline=1, /silent owavl[i] = dum1 nwav[i] = n_elements(dum1) kappa = fltarr(nwav[i], ndust) endelse for iq=0, n_elements(qval_fname)-1 do begin if qval_fname[iq] ne 'mystery' then begin readcol, qval_path[0] + qval_fname[iq], dum1, dum2, skipline=1, /silent readcol, qval_path[0] + qval_fname[iq], nw, gsize, gdens, numline=1, /silent replace = where(dum2 lt 0) if n_elements(replace) gt 0 and replace[0] ne -1 then dum2[replace] = 0.0 ; Interpolate the Qext to the wavelength grid of the data if do_interpolate then begin idum2 = interpol(dum2, dum1, owavl[i]) endif else begin idum2 = interpol(dum2, dum1, owavl[i]) endelse ; Compute mass absorption coefficients from the Qext kappa(*,iq) = 3.0 * idum2 / (4. * gsize[0] * 1d-4 * gdens[0]) print, qval_fname[iq], ' .. reading .... DONE', format='(A40, A30)' endif else begin print, 'Mystery dust detected', format='(A40)' myst_ix = iq endelse endfor kappa_lst.add,kappa endfor ; ; Read the PAH-templates ; ; if npah gt 0 then begin ; pah_prof = fltarr(nwav, npah) ; for ip=0, n_elements(pah_fname)-1 do begin ; readcol, pah_path[0] + pah_fname[ip], dum1, dum2, skipline=1, /silent, format='(F,F)' ; replace = where(dum2 lt 0) ; if n_elements(replace) gt 0 and replace[0] ne -1 then dum2[replace] = 0.0 ; ; Interpolate the PAH profiles to the wavelength grid of the data ; pah_prof[*,ip] = interpol(dum2, dum1, owav) ; print, pah_fname[ip], ' .. reading .... DONE', format='(A40, A30)' ; endfor ; endif return, ndust end function read_all, source_list_file,nr,fit_inc common param_bounds, pbound ;common opacities, kappa_lst, ndust, qval_fname, pah_fname, npah, pah_prof,myst_ix common obs_params, dist, rstar, tstar, nwav, nwav_all, owavl, oflux, ouflux, fit_range, source_name, bll, blpa, rsub, save_chain,lun,fit_incl common grid_params, r_disk_max, nr_disk, geometry, fit_individual_sizes, n_zones common directories, model_dir fit_incl = fit_inc model_dir = './BestFit/' spectra_dir = '' ;'./SPECTRA/' thisname = 'read_all' @natconst.pro ; ; Read the spectral name list ; if fit_incl eq 1 then begin readcol, source_list_file, dum1, dum2, dum3, dum4, dum5, dum6, dum7, dum8, dum9, dum10, dum11,dum12, $ format='(I,A,A,A,F,F,F,F,F,F,F,F)' endif else begin readcol, source_list_file, dum1, dum2, dum3, dum4, dum5, dum6, dum7, dum8, dum9, dum10, dum11, $ format='(I,A,A,A,F,F,F,F,F,F,F)' endelse bll = list() ;list of baseline lengths if fit_incl eq 1 then blpa = list() ;list of baseline position angles owavl = list() ofluxl = list() oufluxl = list() nlines = FILE_LINES(source_list_file) for i = 1,nlines-2 do begin ;1 do begin; ii = where(dum1 eq i) if ii[0] ge 0 then begin if n_elements(ii) gt 1 then begin print, 'ERROR! Multiple files have been found with the same Nr' stop endif else begin spec_fname = dum2(ii(0)) source_name = dum3(ii(0)) sed_fname = dum4(ii(0)) fit_range = [dum5(ii(0)), dum6(ii(0))] rstar = dum7(ii(0)) * rs tstar = dum8(ii(0)) dist = dum9(ii(0)) * pc rsub = dum10(ii(0)) bl = dum11(ii(0)) if fit_incl eq 1 then bla = dum12(ii(0))*!dpi/180.0 endelse endif else begin print, 'ERROR! No spectrum with the specified Nr has been found! ' stop endelse ; ; Read the spectrum itself, select the interval [fit_range[0], fit_range[1]] ; if spec_fname ne 'nopath' then begin readcol, spectra_dir+spec_fname, wav, flux, uflux, skipline=1 ii = where(wav ge fit_range[0] and wav le fit_range[1]) nwav = n_elements(ii) if nwav le 1 then begin print, thisname, ': Not enough datapoints in fitrange' stop endif endif else begin wav = findgen(20,start=fit_range[0],increment=(fit_range[1]-fit_range[0])/(20.0-1.0)) flux = 0.0*wav uflux = 0.0*wav ii = where(wav ge fit_range[0] and wav le fit_range[1]) endelse bll.add,bl if fit_incl eq 1 then blpa.add,bla owavl.add,wav(ii) ofluxl.add,flux(ii) ;print,uflux(ii)/flux(ii)*100.0 ; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ; !!! use nominal 7% error, instead of original data errors! (before it was 5%) oufluxl.add,flux(ii)*0.07 ;uflux(ii) ;flux(ii)*0.05 ;uflux(ii) !!!!!!!!!!!!!!!! endfor nwav_all = 0 nwav = intarr(n_elements(bll)) for i = 0,n_elements(bll)-1 do begin nwav_all = nwav_all + n_elements(owavl[i]) nwav[i] = n_elements(owavl[i]) endfor owav = dblarr(nwav_all) oflux = dblarr(nwav_all) ouflux = dblarr(nwav_all) iw = 0 for i = 0,n_elements(bll)-1 do begin ;print,iw,nwav[i],n_elements(owav[iw:(iw+nwav[i]-1)]) ;print,n_elements(owavl[i]) owav[iw:(iw+nwav[i]-1)] = owavl[i] oflux[iw:(iw+nwav[i]-1)] = ofluxl[i] ouflux[iw:(iw+nwav[i]-1)] = oufluxl[i] iw = iw + nwav[i] endfor return,0 end ; ////////////////////////////////////////////////////////////////////////////////////////// ; Write the results into a file ; ////////////////////////////////////////////////////////////////////////////////////////// pro write_output, ci, comps, spar, redchi, logfile_name common opacities, kappa_lst, ndust, qval_fname, pah_fname, npah, pah_prof,myst_ix common obs_params, dist, rstar, tstar, nwav, nwav_all, owavl, oflux, ouflux, fit_range, source_name, bll, blpa, rsub, save_chain,lun,fit_incl common grid_params, r_disk_max, nr_disk, geometry, fit_individual_sizes, n_zones @natconst ; openw, 1, 'fitresults.log', /append openw, 1, logfile_name, /append norm1 = total(ci(2:ndust+1)) norm2 = total(ci(ndust+2:2*ndust+1)) if n_zones eq 3 then norm3 = total(ci(2*ndust+2:3*ndust+1)) sti = 4 + n_zones if fit_individual_sizes eq 1 then begin q_surfdens1 = spar[sti:(sti+ndust-1)] q_surfdens2 = spar[(sti+ndust):(sti+2*ndust-1)] if n_zones eq 3 then q_surfdens3 = spar[(sti+2*ndust):(sti+3*ndust-1)] endif else begin q_surfdens1 = dblarr(ndust) + spar[sti] q_surfdens2 = dblarr(ndust) +spar[sti+1] if n_zones eq 3 then q_surfdens3 = dblarr(ndust) +spar[sti+2] endelse printf, 1 printf, 1, ' *********************************************** ' printf, 1, ' Input parameters : ' printf, 1, ' *********************************************** ' printf, 1, ' Name of the star : ', source_name printf, 1, ' Radius of the star : ', Rstar/RS printf, 1, ' Temperature of the star : ', Tstar printf, 1, ' Distance : ', dist/pc printf, 1, ' Fitted wavelength range : ', fit_range(0), ' ',fit_range(1) ;printf, 1, ' Sublimation radius : ', r_rim_min/AU ;printf, 1, ' Radius of the rim : ', r_rim_max/AU ;printf, 1, ' Inner disk radius : ', r_disk_min/AU printf, 1, ' Outer disk radius : ', r_disk_max/AU printf, 1, ' *********************************************** ' printf, 1, ' The fitted parameters : ' printf, 1, ' *********************************************** ' ;printf, 1, ' Temperature of the rim : ', spar[0] printf, 1, ' Max temperature of the dust : ', spar[0] ;printf, 1, ' Temperature of the surface l: ', spar[6] printf, 1, ' Max temperature or the cont : ', spar[1] ;printf, 1, ' Prim : ', spar[3] printf, 1, ' pdust : ', spar[2] printf, 1, ' pcont : ', spar[3] ;printf, 1, ' r_gap_in : ', spar[6] printf, 1, ' R1 : ', spar[4] printf, 1, ' R2 : ', spar[5] if n_zones eq 3 then printf, 1, ' R3 : ', spar[6] if fit_incl eq 1 then begin printf, 1, ' cos i : ', spar[-2] printf, 1, ' PA : ', spar[-1]*180.0/!dpi endif printf, 1, ' *********************************************** ' printf, 1, ' Reduced chisquare : ', redchi printf, 1, ' *********************************************** ' printf, 1, ' Contribution of the Star : ', ci[0] ;printf, 1, ' Contribution of the Rim : ', ci[1] printf, 1, ' Contribution of the cont : ', ci[1] printf, 1, ' Contribution of the dust z1 : ', norm1 printf, 1, ' Contribution of the dust z2 : ', norm2 if n_zones eq 3 then printf, 1, ' Contribution of the dust z3 : ', norm3 printf, 1, ' *********************************************** ' printf, 1, ' Dust mass fractions (zone1):' printf, 1, ' *********************************************** ' for i=0, ndust-1 do begin if norm1 eq 0.0 then mf = 0.0 $ else mf = (ci(i+2)/norm1*100.) ;round(ci(i+2)/norm1*100.) printf, 1, qval_fname(i), mf, format='(A40," ", F8.4)' endfor ;if fit_individual_sizes eq 1 then begin printf, 1, ' *********************************************** printf, 1, ' Dust surface density gradients (zone1):' printf, 1, ' *********************************************** ' for i=0, ndust-1 do begin printf, 1, qval_fname(i),q_surfdens1[i], format='(A40," ", F7.3)' endfor ;endif printf, 1, ' *********************************************** ' printf, 1, ' Dust mass fractions (zone2):' printf, 1, ' *********************************************** ' for i=0, ndust-1 do begin if norm2 eq 0. then mf = 0.0 $ else mf = (ci(i+2+ndust)/norm2*100.) ;round(ci(i+2+ndust)/norm2*100.) printf, 1, qval_fname(i), mf, format='(A40," ", F8.4)' endfor ;if fit_individual_sizes eq 1 then begin printf, 1, ' *********************************************** printf, 1, ' Dust surface density gradients (zone2):' printf, 1, ' *********************************************** ' for i=0, ndust-1 do begin printf, 1, qval_fname(i),q_surfdens2[i], format='(A40," ", F7.3)' endfor if n_zones eq 3 then begin printf, 1, ' *********************************************** ' printf, 1, ' Dust mass fractions (zone3):' printf, 1, ' *********************************************** ' for i=0, ndust-1 do begin if norm3 eq 0. then mf = 0.0 $ else mf = (ci(i+2+2*ndust)/norm3*100.) ;round(ci(i+2+2*ndust)/norm3*100.) printf, 1, qval_fname(i), mf, format='(A40," ", F8.4)' endfor ;if fit_individual_sizes eq 1 then begin printf, 1, ' *********************************************** printf, 1, ' Dust surface density gradients (zone3):' printf, 1, ' *********************************************** ' for i=0, ndust-1 do begin printf, 1, qval_fname(i),q_surfdens3[i], format='(A40," ", F7.3)' endfor endif ;endif printf, 1, ' *********************************************** ' printf, 1, ' PAH:' printf, 1, ' *********************************************** ' if npah gt 0 then begin for ip=0, npah-1 do begin printf, 1, pah_fname(ip), ci(ip+3+2*ndust), format='(A40," ", F)' endfor endif close, 1 end ;============================================================================== ; Write an ascii table of the measured spectrum and the various components of ; the best fit model. The columns in the table are as follows: wavelength [mu], ; observed flux [Jy], total model flux [Jy], continuum in model [Jy], PAH [Jy], ; dust [Jy]. Two lines of header, commented with '#' will precede the table. In ; the first line will be the source name, while the second is for the column ; headers. ; Input: model, total model flux ; dust, dust emission from surface layers ; pah, PAH emission ; cont, summed continuum (i.e. midplane, star, ...) ; Return: 0 on success, 1 on error ;------------------------------------------------------------------------------ ; function writeResults, model, dust, pah, cont function writeResults, model, dust, cont common obs_params, dist, rstar, tstar, nwav, nwav_all, owavl, oflux, ouflux, fit_range, source_name, bll, blpa, rsub, save_chain,lun,fit_incl COMMON directories, model_dir fname = model_dir + '/' + source_name + '_results.dat' ; chead = 'wav [mu] flux [Jy] model [Jy] cont. [Jy] PAH [Jy] dust [Jy]' ; format = "(F10.5, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4)" chead = 'wav [mu] flux [Jy] model [Jy] cont. [Jy] dust [Jy]' format = "(F10.5, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4)" openw, fd, fname, /get_lun, width=200 printf, fd, '# ' + source_name printf, fd, '# ' + chead for i = 0, nwav - 1 do begin ; printf,fd, owav[i], oflux[i], model[i], cont[i], pah[i], dust[i], format=format printf,fd, owav[i], oflux[i], model[i], cont[i], dust[i], format=format endfor close, fd free_lun, fd return, 0 end ; ////////////////////////////////////////////////////////////////////////////////////////// ; Planck function ; ////////////////////////////////////////////////////////////////////////////////////////// function planck, wav, temp cc = 2.9979d10 kk = 1.3807d-16 hh = 6.6262d-27 ss = 5.6703d-5 au = 1.496d13 nu = 2.9979d14/wav if temp ne 0.0 then begin I_nu = (2.d*hh*nu^3.d/cc^2) / (exp(hh*nu/(kk*temp))-1) endif else begin if n_elements(wav) > 1 then begin I_nu = dblarr(n_elements(wav)) endif else begin I_nu = 0.0 endelse endelse ;print,hh,cc,kk,temp ;print,'w',wav,'n',nu ;print,'p',planck ;print,'p',temp,min(planck),max(planck),mean(planck) return, I_nu end ; PA [rad] function deproject_uvcoords, pbl, pbla_rad, PA, cos_i ;pbl = sqrt(ucoords^2+vcoords^2) ;pbla_rad = atan(ucoords,vcoords) atd = atan(sin(pbla_rad-PA),(cos(pbla_rad-PA))) ucoords_eff = pbl*(cos(atd)*cos(PA)- cos_i*sin(atd)*sin(PA)) vcoords_eff = pbl*(cos(atd)*sin(PA)+ cos_i*sin(atd)*cos(PA)) pbl_eff = sqrt(ucoords_eff^2+vcoords_eff^2) return, pbl_eff end ; ////////////////////////////////////////////////////////////////////////////////////////// ; Calculate the source function (integrating over radius) ; ////////////////////////////////////////////////////////////////////////////////////////// function get_source_func_rad, wav=wav, spar=spar,do_plot=do_plot common obs_params, dist, rstar, tstar, nwav, nwav_all, owavl, oflux, ouflux, fit_range, source_name, bll, blpa, rsub, save_chain,lun,fit_incl common grid_params, r_disk_max, nr_disk, geometry, fit_individual_sizes, n_zones nw = n_elements(wav) @natconst.pro ;trim0 = spar(0) Tdust0_in = spar[0] Tcont0_in = spar[1] ;prim = spar(3) pdust = spar[2] pcont = spar[3] ;r_gap_in = spar(6)*1.496d13 R1 = spar[4]*1.496d13 ;R2 = spar(5)*1.496d13 ;if n_zones eq 3 then R3 = spar(6)*1.496d13 conv = 1d0 / dist^2. * 1d23 ; ; Create the spatial grid and the corresponding temperatures ; ;r_rim = r_rim_min * (r_rim_max/r_rim_min)^(dindgen(nr_rim)/(nr_rim-1.)) r_disk = R1 * (r_disk_max/R1)^(dindgen(nr_disk)/(nr_disk-1.)) ;print,r_disk/1.496d13 r_disk_rad = r_disk/dist ;r_rim_rad = r_rim/dist ;print,'r_rim',r_rim,'r_disk',r_disk,'g',r_gap_in,R2 ; ; Get the boundary between the out and in temp surf component ; ;rbound = r_disk_min * (Tdust0_out / Tdust0_in)^(1./pdust(0)) ;rbound = r_disk_max ;JV 2022!!! ; ; Compute the temperature grid(s) ; ;if n_zones eq 2 then begin ; idx1 = where(r_disk le R2) ;idx2 = where(r_disk ge R2) ;nr_disk2 = n_elements(idx2) ;nr_disk1 = n_elements(idx1) ;loc1 = value_locate(r_disk, R2) ;!!!!!!!!!!!!! ;loc2 = value_locate(r_disk, R2) ;endif ;if n_zones eq 3 then begin ; idx1 = where(r_disk le R2) ; idx2 = where(r_disk ge R2 and (r_disk lt R3)) ; idx3 = where(r_disk ge R3) ;nr_disk2 = n_elements(idx2) ;nr_disk1 = n_elements(idx1) ;nr_disk3 = n_elements(idx3) ;endif ;print,'nr ',nr_disk,' l ',nr_disk_out,' h ',nr_disk_in,loc1,loc2 ;trim = trim0[0] * (r_rim/r_rim_min)^prim[0] ;t_idx = where(trim lt 3.0,/null) ; ;trim[t_idx] = 0.0 Tdust = Tdust0_in * (r_disk/R1)^pdust ;in ;t_idx = where(Tdust_in lt 3.0,/null) ; ;Tdust_in[t_idx] = 0.0 ;Tdust_out = Tdust_in ;out ;Tdust_out = Tdust0_out * (r_disk/r_disk_min)^pdust ;JV 2019!!! ;Tcont = Tcont0_in * (r_disk/R1)^pcont ;no continuum in this model ;print,'e',Tcont ;t_idx = where(Tcont lt 3.0,/null) ; ;Tcont[t_idx] = 0.0 ; ;plot temperature profile ;print,'u',Tcont ; set_plot,'X' ; window,1 ; r_disk_au = r_disk_rad*180.0/!dpi*3600.0*dist/3.0857200e+18 ; cgplot,r_disk_au,Tdust,thick=th,color=0,xtitle='r (au)',ytitle='T (K)',/ylog ; ; Calculate the emission from the star ; starsurf = !dpi * rstar^2. sflux = planck(wav, tstar) * starsurf * conv ; ; Calculate the source function for each annulus ; ;brim = dblarr(nr_rim,nw) ;brim_arr = dblarr(nr_rim,nw) bcont = dblarr(nr_disk,nw) bcont_arr = dblarr(nr_disk,nw) bdust = dblarr(nr_disk,nw) bdust_arr = dblarr(nr_disk,nw) ; Inner rim ; for ir=0, nr_rim-2 do begin ; da = !dpi * (r_rim(ir+1)*r_rim(ir+1) - r_rim(ir)*r_rim(ir)) ; temp = (trim(ir) + trim(ir+1))*0.5 ; brim_arr(ir+1, *) = planck(wav, temp(0)) * conv ; brim(ir+1, *) = brim_arr(ir+1, *) * da ; endfor ;dust temperature profile ;if(nr_disk_in gt 1) then begin for ir=0, nr_disk-2 do begin da = !dpi * (r_disk[ir+1]^2. - r_disk[ir]^2.) temp = (Tdust[ir] + Tdust[ir+1])*0.5 bdust_arr[ir+1, *] = planck(wav, temp[0]) * conv bdust[ir+1, *] = bdust_arr[ir+1, *]* da endfor ;endif ; out temperature atmosphere ;JV 2019!!! ; if (nr_disk gt nr_disk_in) then begin ; for ir=loc2, nr_disk-2 do begin ; da = !dpi * (r_disk(ir+1)^2. - r_disk(ir)^2.) ; temp = (Tdust_out(ir) + Tdust_out(ir+1))*0.5 ; ;print,ir-loc2-1,temp,(planck(wav, temp[0]) * conv)[10] ; bdust_arr[ir-loc2, *] = planck(wav, temp[0]) * conv ; bdust(ir-loc2, *) = bdust_arr[ir-loc2, *] * da ; endfor ; endif ; continuum temperature profile - no continuum in this model ; for ir=0, nr_disk-2 do begin ; da = !dpi * (r_disk(ir+1)^2. - r_disk(ir)^2.) ; temp = (Tcont(ir) + Tcont(ir+1))*0.5 ; bcont_arr(ir+1, *) = planck(wav, temp(0)) * conv ; bcont(ir+1, *) = bcont_arr(ir+1, *) * da ; endfor ;gap_idx = where((r_disk gt r_gap_in) and (r_disk le R2)) ;bcont_arr[gap_idx,*]=0.0 ;bcont[gap_idx,*]=0.0 ;bcont_arr = bcont ; ; Sum up the contribution of each annulus to each wavelength ; for inu=0, nw-1 do begin ; The rim ;brim(0,inu) = total(brim(1:nr_rim-1,inu), 1) buffer = size(bdust) ; The in temperature atmosphere bdust(0,inu) = total(bdust(1:nr_disk-1,inu), 1) ; The out temperature atmosphere ;if (nr_disk - loc2 gt 1) then bdust(0,inu) = total(bdust(1:*,inu), 1) $ ;JV 2019!!! ;else bdust[0, inu] = 0.d ; The midplane ;bcont(0,inu) = total(bcont(1:nr_disk-1,inu), 1) endfor ;plot if do_plot eq 1 then begin set_plot, 'ps' !P.FONT = 0 r = [0,1,0,0, 1.0,1.0,0.5,0.0,0.5,0.0, 0.0,0.5,0.5, 0.75,0.25,0.25, 0.75,0.50,0.50, 0.75,0.25,0.25,0] * 255 g = [0,0,1,0, 0.5,0.0,1.0,1.0,0.0,0.5, 0.5,0.0,0.5, 0.25,0.75,0.25, 0.50,0.75,0.25, 0.25,0.75,0.50,0] * 255 b = [0,0,0,1, 0.0,0.5,0.0,0.5,1.0,1.0, 0.5,0.5,0.0, 0.25,0.25,0.75, 0.25,0.25,0.75, 0.50,0.50,0.75,0] * 255 tvlct, r,g,b th = 3 plotname = './results_plot/' + source_name + '_source_function.eps' device, filename=plotname, ysize=10.0, xsize=24.0, /color, /encapsulated !p.multi=[0,2,1] plot,wav,bcont[0,*],thick=th,color=0,/ylog oplot,wav,bdust[0,*],thick=th,color=1 ;oplot,wav,bdust[0,*],thick=th,color=2 oplot,wav,sflux,thick=th,color=2 ;oplot,wav,brim[0,*],thick=th,color=4 al_legend, ['bcont', 'bdust','bstar'], linestyle=[0,0,0],$ thick=0*[0,0,0]+th, color=[0,1,2],linsize=0.3,/right_legend idx1=-10 r_disk_au = r_disk_rad*180.0/!dpi*3600.0*dist/3.0857200e+18 plot,r_disk_au,bcont_arr[*,idx1],thick=th,color=0,/ylog,/xlog ;oplot,r_disk_au[loc2:*],bdust_arr[*,idx1],thick=th,color=1 oplot,r_disk_au,bdust_arr,thick=th,color=1 al_legend, ['bcont', 'bdust'], linestyle=[0,0],$ thick=0*[0,0]+th, color=[0,1],linsize=0.3,/right_legend !p.multi=0 device, /close set_plot, 'X' cmd = 'eps2png -f -r 90 '+plotname print,cmd spawn,cmd,cmd_output endif ; ; Create a structure to return ; res = {wav:wav, bstar:sflux, bdust:bdust, bcont:bcont, r_disk_rad:r_disk_rad,$ bdust_arr:bdust_arr, bcont_arr:bcont_arr} return, res end ; ////////////////////////////////////////////////////////////////////////////////////////// ; Scale the paramters from [0,1] to [par_min, par_max] ; ////////////////////////////////////////////////////////////////////////////////////////// function scale_params, par common param_bounds, pbound spar = fltarr(n_elements(par)) ;for i=0, n_elements(par)-1 do begin for i=0, n_elements(pbound(0,*))-1 do begin spar(i) = pbound(0,i) + (pbound(1,i) - pbound(0,i)) * par(i) ; spar[i] = (par[i] - pbound[0,i])/(pbound[1,i] - pbound[0,i]) endfor return, spar end ; ////////////////////////////////////////////////////////////////////////////////////////// ; Scale the paramters from [par_min, par_max] to [0,1] ; ////////////////////////////////////////////////////////////////////////////////////////// function norm_params, spar common param_bounds, pbound par = fltarr(n_elements(spar)) ;print,par ;print,pbound ;for i=0, n_elements(par)-1 do begin for i=0, n_elements(pbound(0,*))-1 do begin par(i)=(spar(i) - pbound(0,i))/(pbound(1,i) - pbound(0,i)) ;print,spar(i),pbound(0,i),pbound(1,i),par(i) if (par[i] ge 1.0) or (par[i] le 0.0) then begin print,'ERROR: parameter ',i,' outside fit range' print,'spar',spar print,'pbound',pbound stop endif endfor return, par end ; ////////////////////////////////////////////////////////////////////////////////////////// ; Calculate the correlated flux (Hankel-transform) ; ////////////////////////////////////////////////////////////////////////////////////////// ; rr [rad] ; pbl_eff [m] ; wl [micron] ; for 1D models function get_corrflux_hankel, radial_profile,rr,pbl_eff,wl common obs_params, dist, rstar, tstar, nwav, nwav_all, owavl, oflux, ouflux, fit_range, source_name, bll, blpa, rsub, save_chain,lun,fit_incl ;rr=rr[-1:*] ;radial_profile=radial_profile[-1:*] ;print,pbl_eff ;b_arg = 2.0*!pi*rr*pbl_eff/(wl*1e-6) ;print,'rr',n_elements(rr),n_elements(radial_profile),min(b_arg),max(b_arg) F_corr = 2.0*!dpi*tsum(rr*dist,dist*rr*radial_profile*beselj(2.0*!dpi*rr*pbl_eff/(wl*1e-6),0)) return, F_corr end ; ////////////////////////////////////////////////////////////////////////////////////////// ; Calculate the model ; ////////////////////////////////////////////////////////////////////////////////////////// function get_model, par, ci=ci, comps=comps,do_bvls=do_bvls,do_plot=do_plot common opacities, kappa_lst, ndust, qval_fname, pah_fname, npah, pah_prof,myst_ix common obs_params, dist, rstar, tstar, nwav, nwav_all, owavl, oflux, ouflux, fit_range, source_name, bll, blpa, rsub, save_chain,lun,fit_incl common grid_params, r_disk_max, nr_disk, geometry, fit_individual_sizes, n_zones common param_bounds, pbound ; ; Scale the parameters ; spar = scale_params(par) ;r_gap_in = spar[6]*1.496d13 pcont = spar[3] R1 = spar[4]*1.496d13 R2 = spar[5]*1.496d13 R1_rad = R1/dist R2_rad = R2/dist if n_zones eq 3 then begin R3 = spar[6]*1.496d13 R3_rad = R3/dist endif ;r_gap_in_rad = r_gap_in/dist if fit_incl eq 1 then begin cos_i = spar[-2] PA = spar[-1] ;print,bll,blpa bll_eff = list() for i=0, n_elements(bll)-1 do begin bll_eff.add,deproject_uvcoords(bll[i], blpa[i], PA, cos_i) endfor endif else begin bll_eff = bll endelse sti = 4 + n_zones if fit_individual_sizes eq 1 then begin q_surfdens1 = spar[sti:(sti+ndust-1)] q_surfdens2 = spar[(sti+ndust):(sti+2*ndust-1)] if n_zones eq 3 then q_surfdens3 = spar[(sti+2*ndust):(sti+3*ndust-1)] endif else begin q_surfdens1 = dblarr(ndust) + spar[sti] q_surfdens2 = dblarr(ndust) +spar[sti+1] if n_zones eq 3 then q_surfdens3 = dblarr(ndust) +spar[sti+2] endelse ;print,'q',q_surfdens,pbound[*,(7+ndust):(7+2*ndust-1)] ; ; Calculate the corresponding source function ; Fdust1 = dblarr(nwav_all, ndust) Fdust2 = dblarr(nwav_all, ndust) if n_zones eq 3 then Fdust3 = dblarr(nwav_all, ndust) bstar = dblarr(nwav_all) ;brim = dblarr(nwav_all) bcont = dblarr(nwav_all) src = get_source_func_rad(wav=owavl[0], spar=spar,do_plot=do_plot) ;print,owavl[0] ; ; Calculate the emissivities ; ;locm = value_locate(src.r_disk_rad, hlr_mid/(dist/3.0857200e+18)/3600.0*!dpi/180.0) ;if locm lt 0 then locm = 0 ;loc1 = value_locate(src.r_disk_rad, R1_rad) loc2 = value_locate(src.r_disk_rad, R2_rad) if n_zones eq 3 then begin loc3 = value_locate(src.r_disk_rad, R3_rad) endif else begin loc3 = -1 endelse ;print,src.r_disk_rad*180.0/!dpi*3600.0*dist/3.0857200e+18 ;print,src.r_disk_rad,'g', r_gap_in_rad,R2_rad ;print,loc2,n_elements(radial_profile2) for idust=0, ndust-1 do begin ;locd = value_locate(src.r_disk_rad, hlr[idust]/(dist/3.0857200e+18)/3600.0*!dpi/180.0) ;if locd lt 0 then locd = 0 ;print,'F_corr',size(F_corr) ; radial_profile1(Nr,Nwl) ; radial_profile1 = ((src.r_disk_rad[0:loc2-1]/src.r_disk_rad[0])^q_surfdens1[idust] # $ ; ((kappa_lst[0])[*,idust])) * src.bdust_arr[0:loc2-1,*] ; radial_profile2 = ((src.r_disk_rad[loc2:(loc3-1)]/src.r_disk_rad[loc2])^q_surfdens2[idust] # $ ; ((kappa_lst[0])[*,idust])) * src.bdust_arr[loc2:(loc3-1),*] ; if n_zones eq 3 then radial_profile3 = ((src.r_disk_rad[loc3:*]/src.r_disk_rad[loc3])^q_surfdens3[idust] # $ ; ((kappa_lst[0])[*,idust])) * src.bdust_arr[loc3:*,*] ; ;print,src.bdust_arr[*,10] ;print,radial_profile2[*,10] ; print,size(src.r_disk_rad) ; print,size(((src.r_disk_rad[0:loc1]/src.r_disk_rad[0])^q_surfdens1[idust] # $ ; ((kappa_lst[i])[*,idust]))) ; print,size(src.bdust_arr) ; print,size(radial_profile1) ; print,size(radial_profile2) ; print,size(src.bdust_arr) iw = 0 ;print,'R12',R1/1.496d13,' ',R2/1.496d13,R3/1.496d13,loc2,loc3,n_elements(src.r_disk_rad),n_elements(radial_profile1[*,0]) ;print,'r',src.r_disk_rad*180.0/!dpi*3600.0*dist/3.0857200e+18 for i=0,n_elements(bll_eff)-1 do begin if myst_ix ne -1 then begin kappa_lst[i,*,myst_ix] = 500.0*(owavl[i]/8.0)^pcont endif radial_profile1 = ((src.r_disk_rad[0:loc2-1]/src.r_disk_rad[0])^q_surfdens1[idust] # $ ((kappa_lst[i])[*,idust])) * src.bdust_arr[0:loc2-1,*] radial_profile2 = ((src.r_disk_rad[loc2:(loc3-1)]/src.r_disk_rad[loc2])^q_surfdens2[idust] # $ ((kappa_lst[i])[*,idust])) * src.bdust_arr[loc2:(loc3-1),*] if n_zones eq 3 then radial_profile3 = ((src.r_disk_rad[loc3:*]/src.r_disk_rad[loc3])^q_surfdens3[idust] # $ ((kappa_lst[i])[*,idust])) * src.bdust_arr[loc3:*,*] ;print,radial_profile ;print,'radial_profile',size(radial_profile) ;print,'bdust_arr',size(src.bdust_arr) ;print,'rr',size(src.r_disk_rad) ;plot,src.r_disk_rad,src.bdust_arr[*,10] ;plot,owavl[i],src.bdust_arr[10,*] ;print,hlr[idust] ;print,'lr',hlr[idust],locd,n_elements(src.r_disk_rad[locd:*]) F_corr1 = dblarr(nwav[i]) F_corr2 = dblarr(nwav[i]) if n_zones eq 3 then F_corr3 = dblarr(nwav[i]) for il=0,nwav[i]-1 do begin ;plot,src.r_disk_rad[0:loc2],radial_profile1[*,il] ;oplot,src.r_disk_rad[loc2:loc3],radial_profile2[*,il],color=255 F_corr1[il] = get_corrflux_hankel(radial_profile1[*,il],src.r_disk_rad[0:loc2-1],bll_eff[i],(owavl[i])[il]) F_corr2[il] = get_corrflux_hankel(radial_profile2[*,il],src.r_disk_rad[loc2:(loc3-1)],bll_eff[i],(owavl[i])[il]) ;------------------------------------------------------------------- ;the next line had an error in summer-autumn 2022: radial_profile2 instead of radial_profile3 if n_zones eq 3 then F_corr3[il] = get_corrflux_hankel(radial_profile3[*,il],src.r_disk_rad[loc3:*],bll_eff[i],(owavl[i])[il]) ;it was radial_profile2 before 2022.11.01 endfor ;print,n_elements(radial_profile1[*,-1]),n_elements(src.r_disk_rad[0:loc2-1]) ;print,owavl[i] ;print,src.r_disk_rad*180.0/!dpi*3600.0*dist/3.0857200e+18 ;plot,owavl[i],F_corr1 ;oplot,owavl[i],F_corr2,color=255 ;plot,owavl[i],radial_profile1[5,*] ;print,'F_corr',size(F_corr) ;print,F_corr ;plot,owavl[i],(kappa_lst[i])[*,idust] * src.bdust(0,*) ;oplot,owavl[i],F_corr,color=255 Fdust1(iw:(iw+nwav[i]-1),idust) = F_corr1 Fdust2(iw:(iw+nwav[i]-1),idust) = F_corr2 if n_zones eq 3 then Fdust3(iw:(iw+nwav[i]-1),idust) = F_corr3 iw = iw+nwav[i] endfor endfor iw = 0 for i=0,n_elements(bll_eff)-1 do begin ;print,(kappa_lst[i])[*,idust] ;print,(kappa_lst[i])[idust,i] * src.bdust(0,*) ;print,Fdust1(iw:(iw+nwav[i]-1),idust) ;print,exp((2.0*hlr[idust]*bll_eff[i]*!pi^2*1000.0/(dist/3.0857200e+18*owavl[i]*648.0))^2/(4.0*alog(2.0))) ;print,owavl[i] ;print,hlr[idust],bll_eff[i],dist/3.0857200e+18 ;Fdust2(*,idust) = kappa(*,idust) * src.bdust(0,*) bstar[iw] = src.bstar ;F_corr1 = dblarr(nwav[i]) ; F_corr_disk = dblarr(nwav[i]) ; for il=0,nwav[i]-1 do begin ; ;F_corr1[il] = get_corrflux_hankel(src.brim_arr[*,il],src.r_rim_rad,bll_eff[i],(owavl[i])[il]) ; ;print,hlr_mid,locm ; F_corr_disk[il] = get_corrflux_hankel((src.bcont_arr[*,il])[*],src.r_disk_rad[*],bll_eff[i],(owavl[i])[il]) ; endfor ;brim[iw] = src.brim[0,*] ;F_corr1 ;iw:(iw+nwav[i]-1) ;bcont[iw] = F_corr_disk iw = iw+nwav[i] endfor ; ; Prepare for the linear decomposition ; ;-------------------- ;factor = [0.54,0.38,0.07,0.03,0.06] ;factor = [0.5, 0.4,0.03,0.06] ;factor= [ 11.00, 44.00, 15.00, 9.00, 2.00, 19.00] ;factor = factor/total(factor) ;Fdust1_new = dblarr(nwav_all)*0.0 ;print,size(Fdust1_new) ;print,size(Fdust1(*,0)) ;for idust=0,ndust-1 do begin ; Fdust1_new = Fdust1_new + factor[idust]*Fdust1(*,idust) ;endfor ;------------------------ ;print,'f',factor ;print,'fn',Fdust1_new bq = dblarr(nwav_all, n_zones*ndust+npah+2) ;orig: dblarr(nwav_all, 2*ndust+npah+3) ;bq = dblarr(nwav_all, 1+npah+3) ;---------- bq(*,0) = bstar ;bq(*,1) = brim ;bq(*,1) = bcont bq(*,2:ndust+1) = Fdust1 ;bq(*,3) = Fdust1_new ;------------- ;print,size(bq(*,3)) ;print,size(Fdust1_new) bq(*,ndust+2:2*ndust+1) = Fdust2 if n_zones eq 3 then bq(*,2*ndust+2:3*ndust+1) = Fdust3 if npah gt 0 then bq(*,n_zones*ndust+2:n_zones*ndust+npah+1) = pah_prof bound = dblarr(2, n_zones*ndust+npah+2) ;orig:dblarr(2, 2*ndust+npah+3) ;bound = dblarr(2, 1+npah+3) ;--------- ;civ = dblarr(1+npah+3) ;------------ bound(0,*) = 0d0 bound(1,*) = 1d4 bound(0,0) = 0.8 ;star bound(1,0) = 1.2 bound(0,1) = 0.0 ;cont. bound(1,1) = 1.0 ;bound(0,2) = 0. ;bound(1,2) = 1. ; ; Do the linear decomposition ; ;print,size(bq) ;print,size(oflux) ;print,size(bound) ;print,size(civ) bk_bq = bq bk_flux = oflux ;print,bk_bq(*,2) if do_bvls eq 1 then begin bvls, bk_bq, bk_flux, bound, ci, iter=iter, ierr=ierr ;bvls, bk_bq, bk_flux, bound, civ, iter=iter, ierr=ierr ;----------- endif ; ; Calculate the model function ; ;--------------- ;bq_new = dblarr(nwav_all, 1*ndust+npah+3) ;bq_new[*,0:2] = bq[*,0:2] ;ci = dblarr( 1*ndust+npah+3) ;ci[0:2]=civ[0:2] ;for idust=0,ndust-1 do begin ; ci[3+idust] = civ[3]*factor[idust] ; bq_new[*,3+idust] = Fdust1[*,idust] ;print,idust,bq_new[*,3+idust] ;print,'f',civ[3]*factor[idust]*Fdust1[*,idust] ;endfor ;bq=bq_new ;---------------- ;print,'ci',ci ; set_plot, 'X' ; idust=0 ; izone=0 ; iw = 0 ; Fcorr=Fdust1[*,idust]*0.0 ; Ftot = Fdust1[*,idust]*0.0 ; for i=0,n_elements(bll_eff)-1 do begin ; Ftot[iw:(iw+nwav[i]-1),idust] = ci[2+0*ndust+idust]*Fdust1[0:nwav[0]-1,idust] ; iw = iw+nwav[i] ; endfor ; ;for idust=0, ndust-1 do begin ; Fcorr = Fcorr+ci[2+0*ndust+idust]*Fdust1[*,idust] ; Fcorr = Fcorr+ci[2+1*ndust+idust]*Fdust2[*,idust] ; Fcorr = Fcorr+ci[2+2*ndust+idust]*Fdust3[*,idust] ; ;endfor ; ; iw = 0 ; i = 0 ; V = Fcorr/Ftot ; window,0 ; cgplot,owavl[i],V[iw:(iw+nwav[i]-1),idust],$ ; xrange=[8.0,13.0],yrange=[min(V)-0.05,max(V)+0.05],thick=3,linestyle=2,xtitle='Wavelenghth (um)',ytitle='V' ; print,'B',bll_eff[i] ; print,'wav',owavl[i] ; print,'V',V[iw:(iw+nwav[i]-1),idust] ; print,'Fcorr',Fcorr[iw:(iw+nwav[i]-1),idust] ; print,'Ftot',Ftot[iw:(iw+nwav[i]-1),idust] ; iw = iw+nwav[i] ; for i=1,n_elements(bll_eff)-1 do begin ; cgplot,owavl[i],V[iw:(iw+nwav[i]-1),idust],/overplot,thick=3,color=0+i*20 ; iw = iw+nwav[i] ; print,'B',bll_eff[i] ; print,'wav',owavl[i] ; ;print,'V',V[iw:(iw+nwav[i]-1),idust] ; ;print,'Fcorr',Fcorr[iw:(iw+nwav[i]-1),idust] ; ;print,'Ftot',Ftot[iw:(iw+nwav[i]-1),idust] ; endfor if do_plot eq 1 then begin ;plot init set_plot, 'ps' !P.FONT = 0 r = [0,1,0,0, 1.0,1.0,0.5,0.0,0.5,0.0, 0.0,0.5,0.5, 0.75,0.25,0.25, 0.75,0.50,0.50, 0.75,0.25,0.25,0] * 255 g = [0,0,1,0, 0.5,0.0,1.0,1.0,0.0,0.5, 0.5,0.0,0.5, 0.25,0.75,0.25, 0.50,0.75,0.25, 0.25,0.75,0.50,0] * 255 b = [0,0,0,1, 0.0,0.5,0.0,0.5,1.0,1.0, 0.5,0.5,0.0, 0.25,0.25,0.75, 0.25,0.25,0.75, 0.50,0.50,0.75,0] * 255 tvlct, r,g,b th = 3 plotname = './results_plot/' + source_name + '_model_plots.eps' device, filename=plotname, ysize=18.0, xsize=24.0,yoffset=24, $ xoffset=0, /color, /landscape, /encapsulated !p.multi=0 first_plot_flag = 0 ;-------------------------------------------------------------------------------------------- ; plot the radial brightness profiles ;-------------------------------------------------------------------------------------------- conv = (1d0 / dist^2. * 1d23) r_disk_au = src.r_disk_rad*180.0/!dpi*3600.0*dist/3.0857200e+18 idx1 = 6 radp_sum1 = 0.0*radial_profile1[*,idx1] radp_sum2 = 0.0*radial_profile2[*,idx1] if n_zones eq 3 then radp_sum3 = 0.0*radial_profile3[*,idx1] for idust=0, ndust-1 do begin ; plot radial surface brightness profiles at a given wavelength ;print,'f',i,idust,il,first_plot_flag ;print,r_disk_au[0:loc1],radial_profile1[*,idx1] radp_sum1 = radp_sum1 + ci[2+idust]*radial_profile1[*,idx1]/conv radp_sum2 = radp_sum2 + ci[2+ndust+idust]*radial_profile2[*,idx1]/conv if n_zones eq 3 then radp_sum3 = radp_sum3 + ci[2+2*ndust+idust]*radial_profile3[*,idx1]/conv endfor if n_zones eq 3 then begin ymin=min([min(radp_sum1),min(radp_sum2),min(radp_sum3)]) ymax=max([max(radp_sum1),max(radp_sum2),max(radp_sum3)]) endif if n_zones eq 2 then begin ymin=min([min(radp_sum1),min(radp_sum2)]) ymax=max([max(radp_sum1),max(radp_sum2)]) endif cgplot,r_disk_au[0:loc2-1],radp_sum1,thick=2*th,color=0,$ ;yrange=[1e-8*max(radp_sum1),max(radp_sum1)],$ /ylog,/xlog,xrange=[0.8*min(r_disk_au),max(r_disk_au)],yrange=[ymin,ymax] cgplot,r_disk_au[loc2:loc3-1],radp_sum2,thick=2*th,color=0,/overplot cgplot,r_disk_au[loc3:*],radp_sum3,thick=2*th,color=0,/overplot for idust=0, ndust-1 do begin cgplot,r_disk_au[0:loc2-1],ci[2+idust]*radial_profile1[*,idx1]/conv,thick=th,color=idust+1,/overplot cgplot,r_disk_au[loc2:loc3-1],ci[2+ndust+idust]*radial_profile2[*,idx1]/conv,thick=th,color=idust+1,linestyle=2,/overplot if n_zones eq 3 then cgplot,r_disk_au[loc3:*],ci[2+2*ndust+idust]*radial_profile3[*,idx1]/conv,thick=th,color=idust+1,linestyle=3,/overplot endfor cgplot,r_disk_au,src.bcont_arr[*,idx1]/conv,thick=th,color=0,linestyle=2,/overplot ;print,src.bcont_arr[*,idx1] dust_labels=[] dust_colors=[] dust_linestyles=[] for j=0,ndust-1 do BEGIN dust_labels=[dust_labels,[qval_fname[j]]] dust_colors=[dust_colors,[1+j]] dust_linestyles = [dust_linestyles,[0]] endfor dust_thick = [2*th,0*dust_colors+th,th] ;print,dust_thick dust_labels=['total dust @'+string((owavl[0])[idx1],format='(%"%6.2f um")'),dust_labels,'cont.'] dust_colors=[0,dust_colors,0] dust_linestyles = [[0],dust_linestyles,[2]] al_legend, dust_labels, linestyle=dust_linestyles,$ thick=dust_thick, color=dust_colors,linsize=0.3,/right_legend !p.multi=0 device, /close set_plot, 'X' cmd = 'eps2png -f -r 90 '+plotname print,cmd spawn,cmd,cmd_output ;-------------------------------------------------------------------------------------------- ; plot the surface density profiles ;-------------------------------------------------------------------------------------------- ; set_plot, 'ps' ; !P.FONT = 0 ; ; r = [0,1,0,0, 1.0,1.0,0.5,0.0,0.5,0.0, 0.0,0.5,0.5, 0.75,0.25,0.25, 0.75,0.50,0.50, 0.75,0.25,0.25,0] * 255 ; g = [0,0,1,0, 0.5,0.0,1.0,1.0,0.0,0.5, 0.5,0.0,0.5, 0.25,0.75,0.25, 0.50,0.75,0.25, 0.25,0.75,0.50,0] * 255 ; b = [0,0,0,1, 0.0,0.5,0.0,0.5,1.0,1.0, 0.5,0.5,0.0, 0.25,0.25,0.75, 0.25,0.25,0.75, 0.50,0.50,0.75,0] * 255 ; tvlct, r,g,b ; ; plotname = './results_plot/' + source_name + '_surfdens_profiles.eps' ; device, filename=plotname, ysize=8.0, xsize=10.0,yoffset=0, $ ; xoffset=0, /color, /portrait, /encapsulated ; ; ;kappa = fltarr(nwav[i], ndust) ; !P.Multi = [0, 1, 1] ; for idust=0,ndust-1 do begin ; sdens_profile1 = src.r_disk_rad[0:loc2-1]/src.r_disk_rad[0])^q_surfdens1[idust] ; sdens_profile2 = src.r_disk_rad[loc2:(loc3-1)]/src.r_disk_rad[loc2])^q_surfdens2[idust] ; if n_zones eq 3 then sdens_profile3 = src.r_disk_rad[loc3:*]/src.r_disk_rad[loc3])^q_surfdens3[idust] ; ; ; endfor ; ; cgplot,r_disk_au[0:loc2-1],radp_sum1,thick=2*th,color=0,yrange=[1e-8*max(radp_sum1),max(radp_sum1)],$ ; /ylog,/xlog,xrange=[0.8*min(r_disk_au),max(r_disk_au)] ; cgplot,r_disk_au[loc2:loc3-1],radp_sum2,thick=2*th,color=0,/overplot ; cgplot,r_disk_au[loc3:*],radp_sum3,thick=2*th,color=0,/overplot ; ; for idust=0, ndust-1 do begin ; cgplot,r_disk_au[0:loc2-1],ci[2+idust]*radial_profile1[*,idx1]/conv,thick=th,color=idust+1,/overplot ; cgplot,r_disk_au[loc2:loc3-1],ci[2+ndust+idust]*radial_profile2[*,idx1]/conv,thick=th,color=idust+1,linestyle=2,/overplot ; if n_zones eq 3 then cgplot,r_disk_au[loc3:*],ci[2+2*ndust+idust]*radial_profile3[*,idx1]/conv,thick=th,color=idust+1,linestyle=3,/overplot ; endfor ; ; ; ; !P.Multi = 0 ; ; ; device, /close ; set_plot, 'X' ; ; cmd = 'eps2png -f '+plotname ; print,cmd ; spawn,cmd,cmd_output ; print,cmd_output endif ;print,size(ci) fmod = bq # ci ; ; Create a structure to return ; comps = bq res = {fmod:fmod, ci:ci} return, res end ; ////////////////////////////////////////////////////////////////////////////////////////// ; Calculate the chisquare ; ////////////////////////////////////////////////////////////////////////////////////////// function get_chi, n, par do_bvls=1 common obs_params, dist, rstar, tstar, nwav, nwav_all, owavl, oflux, ouflux, fit_range, source_name, bll, blpa, rsub, save_chain,lun,fit_incl common opacities, kappa, ndust, qval_fname, pah_fname, npah res = get_model(par,do_bvls=do_bvls,do_plot=0) ;print,oflux ;print,res.fmod ;print,ouflux chi = total((oflux - res.fmod)^2. / ouflux^2.) ; We fit 2*ndust + npah + 3 parameters ci and the BB temperatures of rim, midplane ; and the two atmosphere components + ndust sizes dof = nwav_all - (2*ndust + npah + 3 + 4+ndust) rchi = chi / dof if save_chain then begin printf,lun,rchi,format='(F10.4,$)' for i=0,n_elements(res.ci)-1 do begin printf,lun,string(9B),res.ci[i],format='(A,F14.8,$)' endfor spar = scale_params(par) for i=0,n_elements(spar)-1 do begin printf,lun,string(9B),spar[i],format='(A,F9.4,$)' endfor printf,lun,'' endif return, 1./rchi end ; ////////////////////////////////////////////////////////////////////////////////////////// ; PIKAIA driver routine ; ////////////////////////////////////////////////////////////////////////////////////////// function specfit_tdga, nr=nr, source_list_file=source_list_file, param_file=param_file, logfile_name=logfile_name, rsubl=rsubl, geom=geom, fit_indiv=fit_indiv, n_zon=n_zon, save_chai=save_chai, lu=lu,fit_inc=fit_inc common param_bounds, pbound common opacities, kappa_lst, ndust, qval_fname, pah_fname, npah, pah_prof,myst_ix common obs_params, dist, rstar, tstar, nwav, nwav_all, owavl, oflux, ouflux, fit_range, source_name, bll, blpa, rsub, save_chain,lun,fit_incl common grid_params, r_disk_max, nr_disk, geometry, fit_individual_sizes, n_zones COMMON directories, model_dir @natconst.pro geometry=geom fit_individual_sizes = fit_indiv n_zones = n_zon save_chain = save_chai lun = lu fit_incl = fit_inc if not keyword_set(nr) then nr=1 ;if not keyword_set(source_list_file) then source_list_file='source_list_fit_fulap.inp' if not keyword_set(logfile_name) then logfile_name = 'fitresults.log' ; ; Initialize the fitting reading everything ; dummy=read_all(source_list_file,nr,fit_inc) ; ; Read the opacity and PAH template file names ; ndust = read_opac('opacity_files.inp',1) ; ; Read the parameter input file ; par = read_param_file(param_file,ndust,geometry,fit_individual_sizes,n_zones) pbound = par.pbound ;rin = par.rin ;rrim = par.rrim ;rin = rsub ;rrim = rsub*2.0 rout = par.rout ;nr_rim = par.nr_rim nr_disk = par.nr_disk nvar = par.nvar ctrl = par.ctrl ; ; Spectrum we want to fit ; ;r_rim_min = rsub * au ;rin * au ;r_rim_max = rrim * au ;r_disk_min = rin * au ;rrim * au r_disk_max = rout * au !EXCEPT=0 ;!!!!!!!!!!! ;-------------------------------------------------------------------------------------------- ; RUN THE GENETIC FITTING ROUTINE TO OBTAIN THE BEST TEMPERATURE DISTRIBUTIONS ;-------------------------------------------------------------------------------------------- run_pikaia=1 if run_pikaia eq 1 then begin pikaia, nvar, 'get_chi', ctrl, bestpar, fitness, status model = get_model(bestpar, ci=ci, comps=comps,do_bvls=1,do_plot=1) spar = scale_params(bestpar) endif else begin ;ci =[ 0.80000001, 0.27217068 , 6.7750130e-33 , 3.1019754e-05,$ ;2.9490270e-05, 2.5506162e-05 , 1.5109039e-06, 4.4091475e-06 , 0.0000000,$ ; 0.0000000 , 0.0000000 , 0.0000000 , 0.0000000] ;spar =[ 1477.10 , 1483.96 , 1093.83 , -2.17304 , -0.604095,$ ;-0.209169, 0.0676040 , 0.407306 , 1.5 , 1.5 , 0.4,$ ; 0.41195 , 0.4] ; ci =[ 1.2000000 , 0.32844612 , 4.1442691e-28 , 0.00011273381,$ ;2.4032246e-05 , 7.1179671e-05 , 1.0761427e-05 ,3.4909367e-05 , 0.0000000,$ ; 0.0000000 , 0.0000000 , 0.0000000 , 0.0000000] ; spar =[ 1452.91 , 1439.03 , 1241.72 , -1.74025 , -0.551088,$ ;-0.661711 , 0.356117 , 3.99457 , 0.840973 , 2.04119 , 0.989701,$ ; 0.423896 , -2.28133 , -0.253816 ,-3.47626 , -3.63101 , -1.73861] ;ci =[ 1.2000000 , 0.0000000 , 0.27766181 , 5.3255291e-07,$ ;2.1863210e-06 , 7.2634867e-07 , 3.9620938e-07 , 1.0923872e-07 , 9.1772562e-07] ;spar =[ 0.507348 , 1160.68 , 1498.32 , -1.69606 , -0.300761,$ ;-1.49679 , 3.0 , 1.2, 5.0,2.5,1.2,1.2,1.2] ; ci =[1.0*[ 0.8 , 0.0 , 2e-28],$ ; 0.3*[4e-6 , 4e-7 , 14e-6 , 8e-7, 20e-7, 20e-6 ] ,$ ; 0.3*[10e-6 , 30e-6 , 10e-6 , 40e-6, 16e-7, 16e-6 ]] ; spar =[ 0.515339 , 1500.00 , 1500.0 , -0.918020 , -0.55,$ ; -0.7 , 10.3 , 1.3, -1.5 , -2.9] ; ci =[ 0.80000001 , 0.0000000 , $ ; 0.00010847177 , 3.3633466e-05 ,0.00024889850, 1.0244616e-05 , 0.00016929851 ,$ ; 1.3849242e-05 , 3.3919393e-05 , 0.0000000, 2.7656462e-06 , 7.9666399e-05 ,$ ; 4.4448025e-05 , 0.0000000 , 0.0000000, 2.8022619e-07 , 0.0000000] ; spar =[ 1402.49, 1462.72 , -0.583615 , -0.874452 , 0.208415 , 1.09744,$ ; 2.57925 , -3.60837 , -0.00106597 ,-0.229342] ; Fe: ; ci 0.80000001 0.0000000 4.7730293e-05 1.3457733e-05 2.2082338e-05 ; 4.9859223e-06 0.00013162472 0.00020470117 2.4895917e-05 7.9444958e-05 ; 1.3632679e-05 2.7363849e-06 1.5936806e-05 0.00012606464 1.6644717e-05 ; 2.4406784e-07 0.0000000 7.5349095e-08 4.4506135e-06 0.00010828603 ; spar 1508.08 1.40143 -0.716122 -0.357996 0.175000 1.23946 ; 2.67756 -1.61629 2.50385 1.02857 ; best fit 3 zones, sil + iron (6 dust species) before debug, do not use!!!!!! ci = [ 0.80000001 , 0.000000 , $ 4.3449456e-05, 1.2450783e-05, 1.9481974e-05, 4.5070994e-06, 0.00011828326, 0.00018706673,$ 2.6190348e-05, 9.0018175e-05, 1.5946018e-05, 2.9794365e-06, 1.7238138e-05, 0.00013877550,$ 1.8231803e-05, 5.0756319e-07, 0.0000000, 1.0435137e-07, 4.6636591e-06, 0.00012316086] ;ci = ci*0.0 ;ci[0] = 0.0 ;0.8 ; contribution of the star ;ci[2] = 0.0*4.3449456e-05 ;Q_Am_Mgolivine_Jae_DHS_f0.7_rv2.0.dat zone1 ;ci[14] = 10.0*1.8231803e-05 ;Q_Am_Mgolivine_Jae_DHS_f0.7_rv2.0.dat zone3 spar =[ 1500.0 , 0.01 , -0.729360 , -0.580580,$ 0.180832 , 1.27020 , 2.77823 , -1.50020, $ 2.40011 , 0.913040] ; ; best fit 3 zones, sil + iron (6 dust species), after debug, good fit ; ci 0.80000001 0.0000000 3.0403108e-05 ; 4.9741415e-06 3.0990913e-05 5.4388673e-06 ; 0.00012707724 0.00017523955 1.9628164e-05 ; 2.0824475e-05 6.1674886e-06 7.7430806e-07 ; 9.7302061e-07 7.0569418e-05 0.0000000 ; 0.00086576159 0.0000000 1.7470830e-05 ; 0.00078960547 0.012026179 ; spar 1472.22 0.432140 -0.745225 -0.733773 ; 0.260701 1.00501 4.27855 -3.04059 ; 1.54576 -6.42232 ; bestpar=norm_params(spar) model = get_model(bestpar, ci=ci, comps=comps,do_bvls=0,do_plot=1) endelse chi = total((oflux - model.fmod)^2. / ouflux^2.) dof = nwav_all - (2*ndust + 3 + 4 + ndust) rchi = chi / dof ;-------------------------------------------------------------------------------------------- ; WRITE THE RESULTS INTO A FILE ;-------------------------------------------------------------------------------------------- print,'ci',ci print,'spar',spar for i=0, n_elements(bll)-1 do begin print,'bll (m) = ',bll[i], ' bll_eff (m) = ',deproject_uvcoords(bll[i], blpa[i],spar[-1], spar[-2]),' blpa (deg) = ',blpa[i] endfor ;if run_pikaia eq 1 then write_output, ci, comps, spar, rchi, logfile_name ;-------------------------------------------------------------------------------------------- ; plot the opacity curves ;-------------------------------------------------------------------------------------------- set_plot, 'ps' !P.FONT = 0 r = [0,1,0,0, 1.0,1.0,0.5,0.0,0.5,0.0, 0.0,0.5,0.5, 0.75,0.25,0.25, 0.75,0.50,0.50, 0.75,0.25,0.25,0] * 255 g = [0,0,1,0, 0.5,0.0,1.0,1.0,0.0,0.5, 0.5,0.0,0.5, 0.25,0.75,0.25, 0.50,0.75,0.25, 0.25,0.75,0.50,0] * 255 b = [0,0,0,1, 0.0,0.5,0.0,0.5,1.0,1.0, 0.5,0.5,0.0, 0.25,0.25,0.75, 0.25,0.25,0.75, 0.50,0.50,0.75,0] * 255 tvlct, r,g,b plotname = './results_plot/' + source_name + '_opacity_curves.eps' device, filename=plotname, ysize=ndust*7.0, xsize=12.0,yoffset=0, $ xoffset=0, /color, /portrait, /encapsulated ;kappa = fltarr(nwav[i], ndust) !P.Multi = [0, 1, ndust] for j=0,ndust-1 do begin ii = 0 cgPlot,owavl[ii],(kappa_lst[ii])[*,j],thick=th,xtitle='Wavelenghth (um)',ytitle='kappa_abs [cm2 g-1]' for ii=1,n_elements(bll)-1 do begin cgPlot,owavl[ii],(kappa_lst[ii])[*,j],/overplot endfor al_legend,[qval_fname[j]] endfor !P.Multi = 0 device, /close set_plot, 'X' cmd = 'eps2png -f '+plotname print,cmd spawn,cmd,cmd_output print,cmd_output ;-------------------------------------------------------------------------------------------- ; PLOT THE RESULTS ;-------------------------------------------------------------------------------------------- iw = 0 for ii=0,n_elements(bll)-1 do begin fcomps = fltarr(nwav[ii],n_zones*ndust+2+npah) for i=0, n_zones*ndust+2+npah-1 do begin ;orig 2*ndust+2+npah fcomps(*,i) = comps(iw:(iw+nwav[ii]-1),i) * ci(i) ;print,'ci',i,ci[i] endfor sflux1 = total(fcomps(*,2:ndust+1), 2) sflux2 = total(fcomps(*,ndust+2:2*ndust+1), 2) if n_zones eq 3 then sflux3 = total(fcomps(*,2*ndust+2:3*ndust+1), 2) cont = total(fcomps(*,0:1), 2) ; pahs = total(fcomps(*,2*ndust+3:2*ndust+2+npah),2) norm1 = total(ci(2:ndust+1)) norm2 = total(ci(ndust+2:2*ndust+1)) mfrac1 = ci(2:ndust+1)/norm1*100. mfrac2 = ci(ndust+2:2*ndust+1)/norm2*100. if n_zones eq 3 then begin norm3 = total(ci(2*ndust+2:3*ndust+1)) mfrac3 = ci(2*ndust+2:3*ndust+1)/norm3*100. endif set_plot, 'ps' !P.FONT = 0 r = [0,1,0,0, 1.0,1.0,0.5,0.0,0.5,0.0, 0.0,0.5,0.5, 0.75,0.25,0.25, 0.75,0.50,0.50, 0.75,0.25,0.25,0] * 255 g = [0,0,1,0, 0.5,0.0,1.0,1.0,0.0,0.5, 0.5,0.0,0.5, 0.25,0.75,0.25, 0.50,0.75,0.25, 0.25,0.75,0.50,0] * 255 b = [0,0,0,1, 0.0,0.5,0.0,0.5,1.0,1.0, 0.5,0.5,0.0, 0.25,0.25,0.75, 0.25,0.25,0.75, 0.50,0.50,0.75,0] * 255 tvlct, r,g,b plotname = './results_plot/' + source_name + string(bll[ii],format='(%"_%dm.eps")') device, filename=plotname, ysize=18.0, xsize=24.0,yoffset=24, $ xoffset=0, /color, /landscape, /encapsulated th=3.0 ymin_o = min(oflux[iw:(iw+nwav[ii]-1)]) ymax_o = max(oflux[iw:(iw+nwav[ii]-1)]) ymin_d = min(fcomps[*,2:(2+n_zones*ndust-1)]) ymax_d = max(fcomps[*,2:(2+n_zones*ndust-1)]) ymin_m = min(model.fmod[iw:(iw+nwav[ii]-1)]) ymax_m = max(model.fmod[iw:(iw+nwav[ii]-1)]) ymin_c = min(cont) ymax_c = max(cont) ymin_s1 = min(sflux1) ymax_s1 = max(sflux1) ymin_s2 = min(sflux2) ymax_s2 = max(sflux2) if n_zones eq 3 then begin ymin_s3 = min(sflux3) ymax_s3 = max(sflux3) ymin_s = min([ymin_s1,ymin_s2,ymin_s3]) ymax_s = max([ymax_s1,ymax_s2,ymax_s3]) endif else begin ymin_s = min([ymin_s1,ymin_s2]) ymax_s = max([ymax_s1,ymax_s2]) endelse ymin = min([ymin_o,ymin_d,ymin_m,ymin_c,ymin_s,0.0]) ymax = max([ymax_o,ymax_d,ymax_m,ymax_c,ymax_s]) plot, owavl[ii], oflux[iw:(iw+nwav[ii]-1)], psym=1, title=source_name,thick=th,yrange=[ymin,ymax] ;it was psym=10 oplot, owavl[ii], model.fmod[iw:(iw+nwav[ii]-1)], color=1,thick=2*th oplot, owavl[ii], cont, color=2,thick=3*th oplot, owavl[ii], sflux1, color=1,thick=3*th,linestyle = 2 oplot, owavl[ii], sflux2, color=1,thick=3*th,linestyle= 3 ;sflux_out[iw:(iw+nwav[ii]-1)] if n_zones eq 3 then oplot, owavl[ii], sflux3, color=1,thick=3*th,linestyle=1 dust_labels =[] dust_colors =[] dust_linestyles =[] for j=0,ndust-1 do BEGIN dust_labels=[dust_labels,[qval_fname[j]]+' z1'] dust_colors=[dust_colors,[3+j]] dust_linestyles=[dust_linestyles,[2]] oplot, owavl[ii], fcomps(*,2+j), color=3+j,thick=th,linestyle=2 endfor for j=0,ndust-1 do BEGIN dust_labels=[dust_labels,[qval_fname[j]]+' z2'] dust_colors=[dust_colors,[3+j]] dust_linestyles=[dust_linestyles,[3]] oplot, owavl[ii], fcomps(*,2+ndust+j), color=3+j,thick=th,linestyle=3 endfor if n_zones eq 3 then begin for j=0,ndust-1 do BEGIN dust_labels=[dust_labels,[qval_fname[j]]+' z3'] dust_colors=[dust_colors,[3+j]] dust_linestyles=[dust_linestyles,[1]] oplot, owavl[ii], fcomps(*,2+2*ndust+j), color=3+j,thick=th,linestyle=1 endfor endif ; oplot, owav, pahs, color=7 ; al_legend, ['model', 'cont.', 's_in', 's_out', 'pah'], linestyle=[0,0,0,0,0], color=[2,3,4,6,7] if n_zones eq 2 then al_legend, [['model', 'cont.', 'dust z1','dust z2'],dust_labels], $ linestyle=[[0,0,2,3],dust_linestyles],$ thick=[[2,3,3,3],0*dust_colors+1]*th, color=[[1,2,1,1],dust_colors],linsize=1.0,/right_legend if n_zones eq 3 then al_legend, [['model', 'cont.', 'dust z1','dust z2','dust z3'],dust_labels], $ linestyle=[[0,0,2,3,1],dust_linestyles],$ thick=[[2,3,3,3,3],0*dust_colors+1]*th, color=[[1,2,1,1,1],dust_colors],linsize=1.0,charsize=0.6,/right_legend device, /close set_plot, 'X' cmd = 'eps2png -f -r 90 '+plotname print,cmd spawn,cmd,cmd_output print,cmd_output ; dump = writeResults(model.fmod, sflux_in + sflux_out, pahs, cont) ;dump = writeResults(model.fmod, sflux_in + sflux_out, cont) ;write results fname = model_dir + '/' + source_name + string(bll[ii],format='(%"_%dm_results.dat")') ; chead = 'wav [mu] flux [Jy] model [Jy] cont. [Jy] PAH [Jy] dust [Jy]' ; format = "(F10.5, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4)" chead = ' wav_mu flux_Jy model_Jy cont_Jy dust_Jy ' format = "(F10.5, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ', E12.4, ' ',$)" openw, fd, fname, /get_lun, width=200 printf, fd, '# ' + source_name printf, fd, '# ' + chead, format="(A,$)" for j=0,ndust-1 do BEGIN printf,fd,qval_fname[j],format="(A,' ',$)" endfor printf,fd,'' for i = 0, nwav[ii] - 1 do begin ; printf,fd, (owavl[ii])[i], oflux[i], model[i], cont[i], pah[i], dust[i], format=format printf,fd, (owavl[ii])[i], (oflux[iw:(iw+nwav[ii]-1)])[i], (model.fmod[iw:(iw+nwav[ii]-1)])[i], cont[i], sflux1[i], format=format for j=0,ndust-1 do BEGIN printf,fd, fcomps(i,2+j+2), format="(E12.4, ' ',$)" endfor printf,fd,'' endfor close, fd free_lun, fd iw=iw+nwav[ii] endfor end