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