function ovsa_fit, freq, flux, fmin, fmax, param = param, $ print = print, inter = inter, yrange = yrange ; Define the arrays to be used within the function. ; Theses will be overwritten if fmin and fmax are defined ; or if the INTER keyword is used. freq_new = freq flux_new = flux ; Allow user to select the frequency range overwhich to fit if ; the inter keyword is set jump: if ( keyword_set( inter ) ) then begin pm = !p.multi ; save users p.multi pw = !d.window ; and window settings window, /free circle_sym, /fill if ( keyword_set( yrange ) eq 0 ) then yrange = [ min( flux ), max( flux ) ] plot_oo, freq, flux, /xs, $ xtitle = 'Frequency/GHz', $ ytitle = 'Total Power/SFU', $ psym = 8, yrange = yrange, xrange = [ 1, 20 ], $ title = 'OVSA Spectrum Fitting - select lower then upper frequency' oplot, freq, flux, psym = -1, symsize = 2 print,' ' print,'OVSA> Use the cursor to select a frequency range to fit, ' print,'OVSA> starting with the lower value.' print,' ' cursor, fmin, dum, /data, /down fimin = max( where( freq le fmin ) ) plots, [ freq( fimin ), freq( fimin ) ], $ [ 10^min( !y.crange ), 10^max( !y.crange ) ], linestyle = 0 cursor, fmax, dum, /data, /down fimax = max( where( freq le fmax ) ) plots, [ freq( fimax ), freq( fimax ) ], $ [ 10^min( !y.crange ), 10^max( !y.crange ) ], linestyle = 0 freq_new = freq( fimin : fimax ) flux_new = flux( fimin : fimax ) endif ; If max and min fequency defined then extract relavent sections ; of freq and flux arrays. if ( n_elements( fmin ) gt 0 ) then begin fimin = max( where( freq le fmin ) ) ; closest lower value fimax = min( where( freq ge fmax ) ) ; closest upper value if ( fimin eq -1 ) then fimin = 0 if ( fimax eq -1 ) then fimax = max( where( freq eq max( freq ) ) ) ; set min and/or max to ; 1.2 or 18 if outside ; OVSA freq range freq_new = freq( fimin : fimax ) flux_new = flux( fimin : fimax ) endif ; General initial guess start = [ 1.0, 1.0, 1.0, 1.0 ] ; Define function of the form A x^alpha (1 - exp(-B (x^-beta))) funct = '(p( 0 ) * x^p( 1 )) * ( 1 - exp( -p( 2 )*x^( -p( 3 ) ) ) )' ; Do the fitting result = mpfitexpr( funct, freq_new, flux_new, err, start , /quiet ) fit = fltarr( 2, n_elements( flux_new ) ) fit( 0, * ) = freq_new fit( 1, * ) = result( 0 ) * freq_new^result( 1 ) * $ ( 1 - exp( -result( 2 ) * freq_new^( -result( 3 ) ) ) ) if ( keyword_set( inter ) ) then begin oplot, fit( 0, * ), fit( 1, *), thick = 2 endif ; If the fit is not good, retry. print,' ' ans = '' read,'OVSA> Refit the spectrum (y/n)? ', ans if ( ans eq 'y') then begin wdelete, !d.window !p.multi = pm wset, pw goto, jump endif else begin wdelete, !d.window !p.multi = pm wset, pw endelse ; Store parameters of fit param = [ arr2str( result( 0 ), /trim ), $ arr2str( result( 2 ), /trim ), $ arr2str( result( 1 ), /trim ), $ arr2str( result( 3 ), /trim ), $ arr2str( result( 1 ) - result( 3 ), /trim ) ] ; Print the fit parameters if keyword 'print' set. if ( keyword_set( print ) ) then begin print, ' ' print, 'OVSA> Function form: I = A nu^alpha (1 - exp(-B nu^-beta)) where' print, 'OVSA> A = ' + param( 0 ) print, 'OVSA> B = ' + param( 1 ) print, 'OVSA> alpha = ' + param( 2 ) print, 'OVSA> beta = ' + param( 3 ) print, 'OVSA> alpha - beta = ' + param( 4 ) print, ' ' endif ; Return the fit array return, fit end