PRO conv_spec,sx,sy,px,py,cx,cy

; note
;      no zero padding applied, however sufficient precision also at end points
;      note also, that response function centered, not wrapped arround, 
;      since referred to central frequency of signal.    


n = n_elements(sx)
if n ne n_elements(px) then print,'error: different array dimensions'
shft = long(n/2)
tty = fltarr(n) + (sy(0)+sy(n-1))/2. 
sy_t = sy - tty
cx = sx
cy = shift(fft(fft(sy_t,-1)*fft(py,-1),1),shft)
cy = ((cx(n-1)-cx(0))/(n-1)) * n * cy         ; renormalise result of convolution
cy = cy + tty

END


PRO convol,sx,sy,cx,cy,vdop=vdop,fwhm=fwhm,resol=resol,vsini=vrot,beta=beta,vmacro=vmacro,message=message,original=original,help=help


if n_params() lt 4 then begin 
  
   print,'PRO convol,sx,sy,cx,cy[,vdop=[Km/s],fwhm=,resol=,vsini=[Km/s],beta=,vmacro=[Km/s],/message,/original]'
   if keyword_set(help) then begin
; written by J. Puls, modified by J. Sundqvist 
; spectrum is convolved as specified with keywords. units in 
; Aangstrom or km/s. 

      print,'input'
      print,';           sx          :       wavelength vector'
      print,';           sy          :       flux vector'
      print,';           cx          :       wavelength vector after convolution'
      print,';           cy          :       flux vector after convolution'
;
      print,'; keywords'
      print,';         vdop          :       doppler width [km/s]'
      print,';         fwhm         :        fwhm of gaussian profile [AA] '
      print,';        resol          :       resolution=lambda0/fwhm'   
      print,';        vsini          :       projected rotational velocity [km/s]'
      print,';         beta          :       constant used for rotation'
      print,';       vmacro          :       radial-tangential macroturbulence'
      print,';                               (v_radial=v_tangential) [km/s]'
      print,';      message          :       some information is printed'
      print,';     original          :       if set, remap finally onto original' 
      print,';                               wavelength grid. OBS:linear interpolation.'
;
      print,'; note:'
      print,';      keywords fwhm, resol, vdop mutually exclusive.'
      print,';      if more than one set, the last is ignored'
      print,';      if vsini set, rotational folding is done first'
   endif 
   return
endif
;maximum resolution
max_resol=2d5
;max_resol=3d4

;default settings  
if not keyword_set(beta)     then beta = 1.5 
if not keyword_set(message)  then message = 0
if not keyword_set(original) then original = 0 
if not keyword_set(vrot)     then vrot = 0.d

;constants
clight    = 299792.5d  
rsmpl_fac = 0.7d


;wavelength vector must be strictly increasing! 
i=sort(sx)-lindgen(n_elements(sx))
if max(i) ne 0 then begin
print,'Wavelength vector must be strictly increasing - sort!'
print,'returning'
return
endif

sxx = sx
syy = sy
n   = n_elements(sxx) 
nn  = n
rdst= (sxx(n-1)-sxx(0))/double(n-1)             ; resampling distance
intp= 0


;check if wavelength vector equidistant

dsx        = abs(sxx(1:*)-sxx(0:n-2))
rdst_min   = min(dsx,max=rdst_max)
meanw      = 0.5d*(sxx(0)+sxx(n-1))
eps        = meanw/max_resol
diff       = abs(rdst_max-rdst_min) 
not_eqdist = diff gt eps

if not_eqdist then begin
     if message then print,'spectrum is not equidistant: rdmin= ',rdst_min,', rdmax= ',rdst_max
     rdst = MAX([(rdst_min * rsmpl_fac),eps]) 
     intp = 1
     nn = (sxx(n-1)-sxx(0))/rdst
 endif          

xpow = alog(nn)/alog(2.)
pow = long(xpow)
if ((pow le 14 and pow gt 2) or intp) and ((xpow-pow) ne 0.) then begin
   pow = pow + 1
   nn_x = long(2.^pow)
   if (float(nn)/float(nn_x) gt rsmpl_fac) and (not not_eqdist) then begin
      nn = nn_x * 2l
   endif else begin
      nn = nn_x
   endelse

   rdst = (sxx(n-1)-sxx(0))/(nn-1)

   sxx = lindgen(nn)*rdst + sx(0)

   syy = interpol(sy,sx,sxx,/spline)
   res=finite(syy)
   bres=where(res eq 0, countb)
   if countb ne 0 then begin

   for i=0,nn-1 do begin
   b=where(sx eq sxx(i),count)
     if count gt 1 then begin
       sym=mean(sy(b))
       syy(i)=sym
     endif
   endfor
   endif

   if message then begin
   print,'resampled grid with resampling distance: ',rdst
   print, nn, format= "('resampled spectrum has ',i6,' points')"
endif

endif

px = (findgen(nn)-long(nn/2))*rdst

;rotation

if keyword_set(vrot) and vrot ge 0.1 then begin 
   lambda0 = (sxx(nn-1)+sxx(0))/2.
   dlam_max = vrot*(lambda0/clight)
   xx = px / dlam_max
   one = where(abs(xx) gt 1)
   xx(one) = 1.
   xx_2 = xx*xx
   unsoeld =  (2.*     sqrt(1-xx_2)/!pi + (1-xx_2)*beta/2.) / (1 + 2*beta/3.)
;  gray = (2.*(1-beta)*sqrt(1-xx_2)/!pi + (1-xx_2)*beta/2.) / (1 -   beta/3.) 
   py = unsoeld
   py = py / dlam_max
   conv_spec,sxx,syy,px,py,cx,cy
   sxx = cx
   syy = cy
endif

;macroturbulence

if keyword_set(vmacro) then begin
    aar=1.
    spi=sqrt(!dpi)
    lambda0  = (sxx(nn-1) + sxx(0))/2.
    mr=(vmacro*lambda0/clight)[0]
    ccr=2.*aar/spi/mr
    pxmr=abs(px)/mr
    py=ccr*[exp(-pxmr^2)+spi*pxmr*(errorf(pxmr)-1.)]
    py=py/int_tabulated(px,py,/double); to be sure the distribution is normalized
    conv_spec,sxx,syy,px,py,cx,cy
    sxx = cx
    syy = cy 
endif   
     
;Gauss profile

if keyword_set(vdop) or keyword_set(fwhm) or keyword_set(resol) then begin

   fac     = 2.*sqrt(alog(2.))
   lambda0 = (sxx(nn-1)+sxx(0))/2.
  
  if keyword_set(fwhm) then begin
    dlam_dop = fwhm/fac
    vdop     = dlam_dop*clight/lambda0

  endif else if keyword_set(resol) then begin
    vdop=clight/(resol*fac)
    dlam_dop = vdop*(lambda0/clight)
 
  endif else begin  
    dlam_dop = vdop*(lambda0/clight)
  endelse

  min_cnr = 5                  
  max_cnr = long(nn*0.8)

   if message then begin 
      print,    vdop,    dlam_dop, format="('vdop=',f7.3,'km/s  <=>',f7.3,'a')"
      print,fac*vdop,fac*dlam_dop, format="('fwhm=',f7.3,'km/s  <=>',f7.3,'a')"
   endif
   cnr = long(4.*dlam_dop/rdst) + 1       
   if max_cnr le cnr then begin
      percent = cnr/float(nn)*100.
      print,'gauss profile is very large:' + string(percent,'(f6.2)') $
                               + '% of spectrum points are forming the profile'
    endif
   if min_cnr ge cnr then begin
   print,'gauss profile is very small: only ' + string(cnr,'(i4)') $  
                                       + ' data points are forming the profile'
   endif
   xx = px / dlam_dop
   six = where(abs(xx) lt 6.)
   py = fltarr(nn)
   py(six) = exp(-(xx(six)^2))/(sqrt(!pi)*dlam_dop)
   conv_spec,sxx,syy,px,py,cx,cy       
   sxx = cx
   syy = cy
endif

if original then begin
   cy = interpol(syy,sxx,sx)
   cx = sx
endif

;Return only real part
cy=real_part(cy)
cx=real_part(cx)       

END

function fconvol,x,y,xout,rv=rv,vdop=vdop,fwhm=fwhm,resol=resol,vsini=vrot,beta=beta,vmacro=vmacro,message=message
  xx=x[sort(x)]
  yy=y[sort(x)]
  yy=yy[uniq(xx)]
  xx=xx[uniq(xx)]
  n=n_elements(xx)
  cc=2.99792458d5
  if keyword_set(rv) then xx=xx*(1.0d - rv/cc)
  if n_elements(xx) gt 50000 then begin
     nps=25000l
     segs=(n/nps)    
     yff=yy*0.0d
     if n mod nps lt 1000 then begin
        nps=24000l
        segs=(n/nps)
     endif
     for ii=0, segs do begin 
        ind=[0l,nps-1l]+ii*nps
        if (ii+1l)*nps gt n then ind[1]=n-1
        cind=[max([0,ind[0]-1000l]), min([ind[1]+1000l,n-1])]
        yff[ind[0] : ind[1]]=fconvol(xx[cind[0] : cind[1]],yy[cind[0] : cind[1]],xx[ind[0] : ind[1]],vdop=vdop,fwhm=fwhm,resol=resol,vsini=vrot,beta=beta,vmacro=vmacro,message=message)
     endfor
     yout=interpol(yff,xx,xout)


     
     return,yout
  endif
  if odd(n_elements(xx)) then begin
     xx=xx[1:*]
     yy=yy[1:*]
  endif
  convol,xx,yy,nx,sol,vdop=vdop,fwhm=fwhm,resol=resol,vsini=vrot,beta=beta,vmacro=vmacro,message=message,/original
  if n_elements(xout) gt 1 then yout=interpol(sol,xx,xout) else yout=sol 
  return,yout
end