function gbroad, dw, s, fwhm
;Smooths a spectrum by convolution with a gaussian of specified fwhm.
; dw   (scalar) wavelength/velocity step of the spectrum to be smoothed
; s    (vector) spectrum to be smoothed
; fwhm (scalar) full width at half maximum, in wavelength/velocity units,
;      of smoothing gaussian.
;Returns a vector containing the gaussian-smoothed spectrum.
;History:
; 23-Feb-07 Kochukhov Adapted from gaussbroad.pro

if n_params() lt 3 then begin
  print,'syntax: snew=gbroad(dw,sold,fwhm)'
  retall
endif

;Warn user if hwhm is negative.
  if fwhm lt 0.0 then $
    message,/info,'Warning! Forcing negative smoothing width to zero.'

;Return input argument if half-width is nonpositive.
  if fwhm le 0.0 then return,s			;true: no broadening

;Number of spectral points and hwhm
  nw = n_elements(s)
  hwhm = 0.5 * fwhm

;Make smoothing gaussian; extend to 4 sigma.
;Note: 4.0 / sqrt(2.0*alog(2.0)) = 3.3972872 and sqrt(alog(2.0))=0.83255461
;  sqrt(alog(2.0)/pi)=0.46971864 (*1.0000632 to correct for >4 sigma wings)
  if(hwhm gt 5*dw*(nw-1)) then return, replicate(total(s)/nw, nw)
  nhalf = long(3.3972872 * hwhm/dw)		;# points in half gaussian
  ng = 2 * nhalf + 1				;# points in gaussian (odd!)
  wg = dW * (findgen(nG) - (ng-1)/2.0)		;wavelength scale of gaussian
  xg = (0.83255461 / hwhm) * wg 		;convenient absisca
  gpro = (0.46974832 * dw / hwhm) * exp(-xg*xg)	;unit area gaussian w/ FWHM

;Pad spectrum ends to minimize impact of Fourier ringing.
  npad = nhalf + 2				;# pad pixels on each end
  spad = [replicate(s(0),npad),s,replicate(s(nw-1),npad)]

;Convolve and trim.
  sout = convol(spad,gpro)			;convolve with gaussian
  sout = sout(npad:npad+nw-1)			;trim to original data/length
  return,sout					;return broadened spectrum.

end

function rtint,mu,inten,deltav,vsini_in,vrt_in,osamp=osamp
;+
;NAME:
;	RTINT
;
;PURPOSE:
;	Produces a flux profile by integrating intensity profiles (sampled
;	  at various mu angles) over the visible stellar surface.
;
;CALLING SEQUENCE:
;	flux = RTINT(mu, inten, deltav, vsini, vrt)
;
;INPUTS:
;	MU: (vector(nmu)) cosine of the angle between the outward normal and
;	  the line of sight for each intensity spectrum in INTEN.
;	INTEN: (array(npts,nmu)) intensity spectra at specified values of MU.
;	DELTAV: (scalar) velocity spacing between adjacent spectrum points
;	  in INTEN (same units as VSINI and VRT).
;	VSINI (scalar) maximum radial velocity, due to solid-body rotation.
;	VRT (scalar) radial-tangential macroturbulence parameter, i.e.
;         sqrt(2) times the standard deviation of a Gaussian distribution
;         of turbulent velocities. The same distribution function describes
;	  the radial motions of one component and the tangential motions of
;	  a second component. Each component covers half the stellar surface.
;         See _The Observation and Analysis of Stellar Photospheres_, Gray.
;
;INPUT KEYWORDS:
;	OSAMP: (scalar) internal oversampling factor for convolutions. By
;	  default convolutions are done using the input points (OSAMP=1),
;	  but when OSAMP is set to higher integer values, the input spectra
;	  are first oversampled by cubic spline interpolation.
;
;OUTPUTS:
;	Function Value: (vector(npts)) Disk integrated flux profile.
;
;RESTRICTIONS:
;	Intensity profiles are weighted by the fraction of the projected
;	  stellar surface they represent, apportioning the area between
;	  adjacent MU points equally. Additional weights (such as those
;	  used in a Gauss-Legendre quadrature) can not meaningfully be
;	  used in this scheme.  About twice as many points are required
;	  with this scheme to achieve the precision of Gauss-Legendre
;	  quadrature.
;	DELTAV, VSINI, and VRT must all be in the same units (e.g. km/s).
;	If specified, OSAMP should be a positive integer.
;
;AUTHOR'S REQUEST:
;	If you use this algorithm in work that you publish, please cite
;	  Valenti & Anderson 1996, PASP, currently in preparation.
;
;MODIFICATION HISTORY:
;	   Feb-88  GM	Created ANA version.
;	13-Oct-92 JAV	Adapted from G. Marcy's ANA routine of the same name.
;	03-Nov-93 JAV	Switched to annular convolution technique.
;	12-Nov-93 JAV	Fixed bug. Intensity components not added when vsini=0.
;	14-Jun-94 JAV	Reformatted for "public" release. Heavily commented.
;			Pass deltav instead of 2.998d5/deltav. Added osamp
;			  keyword. Added rebinning logic at end of routine.
;			Changed default osamp from 3 to 1.
;	20-Feb-95 JAV	Added mu as an argument to handle arbitrary mu sampling
;			  and remove ambiguity in intensity profile ordering.
;			Interpret VTURB as sqrt(2)*sigma instead of just sigma.
;			Replaced call_external with call to spl_{init|interp}.
;	03-Apr-95 JAV	Multiply flux by !pi to give observed flux.
;	24-Oct-95 JAV	Force "nmk" padding to be at least 3 pixels.
;	18-Dec-95 JAV	Renamed from dskint() to rtint(). No longer make local
;			  copy of intensities. Use radial-tangential instead
;			  of isotropic Gaussian macroturbulence.
;	26-Jan-99 JAV	For NMU=1 and VSINI=0, assume resolved solar surface;
;			  apply R-T macro, but supress vsini broadening.
;	01-Apr-99 GMH	Use annuli weights, rather than assuming equal area.
;-

if n_params() lt 5 then begin
  print,'syntax: flux = rtint(mu,inten,deltav,vsini,vrt [,osamp=])'
  retall
endif

;Make local copies of various input variables, which will be altered below.
  vsini = float(vsini_in)			;ensure real number
  vrt = float(vrt_in)				;ensure real number

;Determine oversampling factor.
  if n_elements(osamp) eq 0 then osamp = 0	;make sure variable is defined
  os = round(osamp > 1)				;force integral value > 0

;Convert input MU to projected radii, R, of annuli for a star of unit radius
;  (which is just sine, rather than cosine, of the angle between the outward
;  normal and the line of sight).
  rmu = sqrt(1.0 - mu^2)			;use simple trig identity

;Sort the projected radii and corresponding intensity spectra into ascending
;  order (i.e. from disk center to the limb), which is equivalent to sorting
;  MU in descending order.
  isort = sort(rmu)				;sorted indicies
  rmu = rmu(isort)				;reorder projected radii
  nmu = n_elements(mu)				;number of radii 
  if(nmu eq 1) then vsini=0.0                   ;ignore vsini if only 1 mu 

;Calculate projected radii for boundaries of disk integration annuli.  The n+1
;  boundaries are selected such that r(i+1) exactly bisects the area between
;  rmu(i) and rmu(i+1). The innermost boundary, r(0) is set to 0 (disk center)
;  and the outermost boundary, r(nmu) is set to 1 (limb).
  if nmu gt 1 or vsini ne 0 then begin		;really want disk integration
    r = sqrt(0.5 * ( rmu(0:nmu-2) ^ 2 $
                   + rmu(1:nmu-1) ^ 2 ) )	;area midpoints between rmu
    r = [0, temporary(r), 1]			;bookend with center and limb

;Calculate integration weights for each disk integration annulus.  The weight
;  is just given by the relative area of each annulus, normalized such that
;  the sum of all weights is unity.  Weights for limb darkening are included
;  explicitly in the intensity profiles, so they aren't needed here.
    wt = r(1:nmu) ^ 2 - r(0:nmu-1) ^ 2		;weights = relative areas
  endif else begin
    wt = [1.0]					;single mu value, full weight
  endelse

;Generate index vectors for input and oversampled points. Note that the
;  oversampled indicies are carefully chosen such that every "os" finely
;  sampled points fit exactly into one input bin. This makes it simple to
;  "integrate" the finely sampled points at the end of the routine.
  vinfo = size(inten)			;variable information block
  npts  = long(vinfo(1))		;# of points
  xpix  = dindgen(npts)			;point indices
  nfine = long(os*npts)			;# of oversampled points
  xfine = (0.5 / os) $
	* (2*dindgen(nfine)-os+1) 	;oversampled points indices

;Loop through annuli, constructing and convolving with rotation kernels.
  dummy = 0L				;init call_ext() return value
  yfine = dblarr(nfine,/nozero)		;init oversampled intensities
  flux = dblarr(nfine)			;init flux vector
  for imu=0,nmu-1 do begin		;loop thru integration annuli

;Use external cubic spline routine (adapted from Numerical Recipes) to make
;  an oversampled version of the intensity profile for the current annulus.
;  IDL (tensed) spline is nice, but *VERY* slow. Note that the spline extends
;  (i.e. extrapolates) a fraction of a point beyond the original endpoints.
    ypix = double(inten(*,isort(imu)))	;extract intensity profile
    if os eq 1 then begin		;true: no oversampling
      yfine = ypix			;just copy (use) original profile
    endif else begin			;else: must oversample
      secder = spl_init(xpix, ypix)	;find spline coefficients
      yfine = spl_interp(xpix, ypix $
	               ,secder, xfine)	;spline onto fine wavelength scale
    endelse

;Construct the convolution kernel which describes the distribution of
;  rotational velocities present in the current annulus. The distribution has
;  been derived analytically for annuli of arbitrary thickness in a rigidly
;  rotating star. The kernel is constructed in two pieces: one piece for
;  radial velocities less than the maximum velocity along the inner edge of
;  the annulus, and one piece for velocities greater than this limit.
    if vsini gt 0 then begin			;true: nontrivial case
      r1 = r(imu)				;inner edge of annulus
      r2 = r(imu+1)				;outer edge of annulus
      dv = deltav / os				;oversampled velocity spacing
      maxv = vsini * r2				;maximum velocity in annulus
      nrk = 2*long(maxv/dv) + 3			;# oversampled kernel point
      v = dv * (dindgen(nrk) - ((nrk-1)/2))	;velocity scale for kernel
      rkern = dblarr(nrk)			;init rotational kernel
      j1 = where(abs(v) lt vsini*r1,n1)		;low velocity points
      if n1 gt 0 then begin
        rkern(j1) $
	  = sqrt((vsini*r2)^2 - v(j1)^2) $
	  - sqrt((vsini*r1)^2 - v(j1)^2)	;generate distribution
      endif
      j2 = where(abs(v) ge vsini*r1 $		;higher velocity points
             and abs(v) le vsini*r2,n2)
      if n2 gt 0 then begin
        rkern(j2) $
	  = sqrt((vsini*r2)^2 - v(j2)^2)	;generate distribution
      endif
      rkern = rkern / total(rkern)		;normalize kernel

;Convolve the intensity profile with the rotational velocity kernel for this
;  annulus. Pad each end of the profile with as many points as are in the
;  convolution kernel. This reduces Fourier ringing. The convolution may also
;  be done with a routine called "externally" from IDL, which efficiently
;  shifts and adds.
      if nrk gt 3 then begin
;       lpad = replicate(yfine(0),nrk)		;padding for the "left" side
;       rpad = replicate(yfine(nfine-1),nrk)	;padding for the "right" side
;       yfine = convol([lpad,yfine,rpad],rkern)	;add the padding and convolve
;       yfine = yfine(nrk:nrk+nfine-1)		;trim away padding
        yfine = convol(yfine, rkern, /edge_truncate)
      endif
    endif					;done with rotation

;Calculate projected sigma for radial and tangential velocity distributions.
    muval = mu(isort(imu))			;current value of mu
    sigma = os * vrt / sqrt(2) / deltav		;standard deviation in points
    sigr = sigma * muval			;reduce by current mu value
    sigt = sigma * sqrt(1.0 - muval^2)		;reduce by sqrt(1-mu^2)

;Figure out how many points to use in macroturbulence kernel.
    nmk = (fix(10*sigma)) < ((nfine-3)/2)	;extend kernel to 10 sigma
    nmk = nmk > 3				;pad with at least 3 pixels

;Construct radial macroturbulence kernel with a sigma of mu*VRT/sqrt(2).
    if sigr gt 0 then begin
      xarg = (dindgen(2*nmk+1)-nmk) / sigr	;exponential argument
      mrkern = exp((-0.5*(xarg^2)) > (-20))	;compute the gaussian
      mrkern = mrkern / total(mrkern)		;normalize the profile
    endif else begin
      mrkern = dblarr(2*nmk+1)			;init with 0d0
      mrkern(nmk) = 1.0				;delta function
    endelse

;Construct tangential kernel with a sigma of sqrt(1-mu^2)*VRT/sqrt(2).
    if sigt gt 0 then begin
      xarg = (dindgen(2*nmk+1)-nmk) / sigt	;exponential argument
      mtkern = exp((-0.5*(xarg^2)) > (-20))	;compute the gaussian
      mtkern = mtkern / total(mtkern)		;normalize the profile
    endif else begin
      mtkern = dblarr(2*nmk+1)			;init with 0d0
      mtkern(nmk) = 1.0				;delta function
    endelse

;Sum the radial and tangential components, weighted by surface area.
    area_r = 0.5				;assume equal areas
    area_t = 0.5				;ar+at must equal 1
    mkern = area_r * mrkern + area_t * mtkern	;add both components

;Convolve the total flux profiles, again padding the spectrum on both ends to
;  protect against Fourier ringing.
    lpad = replicate(yfine(0),nmk)		;padding for the "left" side
    rpad = replicate(yfine(nfine-1),nmk)	;padding for the "right" side
    yfine = convol([lpad,yfine,rpad],mkern)	;add the padding and convolve
    yfine = yfine(nmk:nmk+nfine-1)		;trim away padding
;   yfine = convol(yfine, mkern, /edge_truncate)

;Add contribution from current annulus to the running total.
    flux = flux + wt(imu)*yfine	                ;add profile to running total

  endfor

  flux = reform(flux,os,npts,/overwrite)	;convert to an array
  return,!pi * total(flux,1) / os		;sum, normalize, and return

end

pro rds3v, infile, outfile, x, spl, cf, $
          vsini=vsini, vrt=vrt, resol=resol, instrum=instrum, $
          nosave=nosave, stokes=stokes, spq=spq, spu=spu, spv=spv, spi=spi, spc=spc

;Inputs:
; infile  (string) name of ".out" file generated by Synth3.
; outfile (string) name of output file in the format of "rotate".
;                  outfile should have no extension, file of
;                  type outfile+'.prf' will be created.
; vsini   (scalar) vsini (in km/s) for disk integration.
; vrt     (scalar) radial-tangential macroturbulence (in km/s).
; resol   (scalar) resolution for spectral rebinning (lambda/delta lambda)
; instrum (scalar) instrumental profile (lambda/delta lambda)
; nosave  (flag)   inhibits writing output file
; stokes  (flag)   integrate polarization spectra
;Outputs:
; x   wavelength scale for output spectrum
; spl flux spectrum constructed from intensities in file.
; cf  formation depth in log tau(5000)
;Requires:
;  rtint.pro
;History:
; 23-Feb-07 Kochukhov Adapted from rdmag.pro

if n_params() lt 2 then begin
  print,'syntax: rds3v, infile, outfile, [w, i, cf, vsini= ,vrt=, resol=, instrum= , /nosave, /stokes]'
  print,'  infile - input Synth3/Synthmag spectrum synthesis file'
  print,' outfile - output file for Rotate/Binmag'
  print,'       w - wavelength array'
  print,'       i - Stokes I array'
  print,'      cf - contribution function array'
  print,'   vsini - projected rotational velocity (km/s)'
  print,'     vrt - radial-tangential macroturbulent velocity (km/s)'
  print,'   resol - velocity sampling of the output spectrum (default lambda/delta_lambda=10^6)'
  print,' instrum - instrumental broadening lambda/delta_lambda'
  print,' /nosave - do not create output spectra'
  print,' /stokes - save output spectra for all Stokes parameters'
  print
  print,' To write Rotate/Binmag spectra: '
  print,'   rds3v,infile,outfile,vsini=vsini,vrt=vrt,instrum=instrum'
  print,' To read in IDL without saving files: '
  print,'   rds3v,infile,'''',w,i,vsini=vsini,vrt=vrt,instrum=instrum,/nosave'
  return
endif

;Set defaults for optional keywords.
  if n_elements(vrt    ) eq 0 then vrt     = 0.0	    ;RT macroturbulence 
  if n_elements(vsini  ) eq 0 then vsini   = 0.0	    ;rotational velocity
  if n_elements(resol  ) eq 0 then resol   = 1d6            ;resolution for output spectrum

;Active nosave if output file array is empty
  if not keyword_set(nosave) and strlen(strtrim(outfile,2)) lt 1 then nosave=1

  goto,open_files

;In case of crash close all files
close_files:
  close,/all
  print,'Error reading input file '+"'"+infile+"'"
  return

;Open files.
;If /nosave then we will not write output to file
open_files:
  ON_IOERROR,close_files
  openr, unit, infile, /get, comp=(strmid(infile,strlen(infile)-2,2) eq 'gz')
  if not keyword_set(nosave) then begin
    openw, outI, outfile+'.prf', /get
    if keyword_set(stokes) then begin
      openw, outQ, outfile+'Q.prf', /get
      openw, outU, outfile+'U.prf', /get
      openw, outV, outfile+'V.prf', /get
    endif
  endif
  
;Copy line list.
  nlines=0l
  readf, unit, nlines 
  if not keyword_set(nosave) then begin
    printf, outI, nlines,' - number of spectral lines'
    if keyword_set(stokes) then begin
      printf, outQ, nlines,' - number of spectral lines'
      printf, outU, nlines,' - number of spectral lines'
      printf, outV, nlines,' - number of spectral lines'
    endif
  endif
  
  sbuff = ''
  for i=1l,nlines do begin
    readf, unit, sbuff
    if not keyword_set(nosave) then begin
      printf, outI, sbuff
      if keyword_set(stokes) then begin
        printf, outQ, sbuff
        printf, outU, sbuff
        printf, outV, sbuff
      endif
    endif  
  endfor

;Predefine variable types for use in read statements below.
  nmu = 0
  nwl = 0L
  nwc = 0
  
;Read mu values  
  readf, unit, nmu
  mu = dblarr(nmu)
  readf, unit, mu

;Read continuum intensities
  readf, unit, nwc
  spc = dblarr(nmu+1,nwc)
  readf, unit, spc
  wlc = reform(spc[0,*])
  spc = spc[1:nmu,*]

;Read line intensities
  readf, unit, nwl
  spl = dblarr(nmu+1,nwl)
  readf, unit, spl
  wll = reform(spl[0,*])
  spl = spl[1:nmu,*]

;Read polarization spectra
  if keyword_set(stokes) then begin
    spq = dblarr(nmu+1,nwl)
    spu = spq
    spv = spq
    readf, unit, spq
    readf, unit, spu
    readf, unit, spv
    spq = spq[1:nmu,*]
    spu = spu[1:nmu,*]
    spv = spv[1:nmu,*]
  endif

;Read formation depth
  if n_params() gt 4 and not eof(unit) then begin
    cf=dblarr(4, nwl)
    readf, unit, cf
  endif
  
;Close input file
  free_lun, unit

  wfirst=wlc[0] & wlast=wlc[nwc-1]
  nx = floor(alog10(wlast/wfirst)/alog10(1d0+1d0/resol))+1
  rstep = 1d0/((wlast/wfirst)^(1d0/(nx-1d0))-1d0)            ; actual resolution
  deltav = 2.99792458d5 / rstep

;Interpolate onto grid equispaced in velocity 
  x = wfirst * (1d0 + 1d0 / rstep)^dindgen(nx)
  sp1l = dblarr(nx,nmu)
  if keyword_set(stokes) then begin
    sp1q = sp1l & sp1u = sp1l & sp1v = sp1l
  endif
  sp1c = dblarr(nx,nmu)
  sp2c = dblarr(nmu)
  for imu = 0, nmu-1 do begin
    sp1l(*,imu) = interpol(spl(imu,*), wll, x)
    if keyword_set(stokes) then begin
      sp1q(*,imu) = interpol(spq(imu,*), wll, x)
      sp1u(*,imu) = interpol(spu(imu,*), wll, x)
      sp1v(*,imu) = interpol(spv(imu,*), wll, x)
    endif
    sp1c(*,imu) = interpol(spc(imu,*), wlc, x)
    sp2c(imu) = interpol(spc(imu,*), wlc, 0.5*(wfirst+wlast))
  endfor  

;Do the disk integration.
  spl = rtint(mu, sp1l, deltav, vsini, vrt)
  if keyword_set(stokes) then begin
    spq = rtint(mu, sp1q, deltav, vsini, vrt)
    spu = rtint(mu, sp1u, deltav, vsini, vrt)
    spv = rtint(mu, sp1v, deltav, vsini, vrt)
  endif
  spc = rtint(mu, sp1c, deltav, vsini, vrt)

;Integrated continuum intensity at the ends of spectral interval
  fcblue = spc[0   ]/2./!dpi
  fcred  = spc[nx-1]/2./!dpi
  
;Apply instrumental broadening
  if n_elements(instrum) gt 0 then begin
    fwhm = 2.99792458d5 / instrum
    spl = gbroad(deltav, spl, fwhm)
    if keyword_set(stokes) then begin
      spq = gbroad(deltav, spq, fwhm)
      spu = gbroad(deltav, spu, fwhm)
      spv = gbroad(deltav, spv, fwhm)
    endif
  endif

;Normalize to continuum
  spi = spl
  spl = spl / spc
  if keyword_set(stokes) then begin
    spq = spq / spc
    spu = spu / spc
    spv = spv / spc
  endif

;Write output in the format of "rotate"
  if not keyword_set(nosave) then begin

;Find intensity in the disk center by qaudratic extrapolation
    c = poly_fit(mu,sp2c,2)
    fc0 = poly([1.0],c)

;Find linear limb darkening coefficient.
    c = poly_fit(mu,sp2c,1)
    limbd = c[1]/poly([1.0],c)
    ;limbd = 1.0 - c[0]/poly([1.0],c)         ; gives identical result

    fcscale=1.0686475D5
    printf, outI, nx, limbd, form='(I10,F8.4," - number of wavelengths, limb darkening")'
    for i=0L,nx-1 do printf, outI, x(i), spl(i)
    printf, outI, ' '+infile
    printf, outI, fcblue*fcscale, fcred*fcscale, fc0*fcscale, $
            form= '(" FCblue=",E12.6," FCred=",E12.6," FC(0)=",E12.6)'
    free_lun, outI
    if keyword_set(stokes) then begin
      printf, outQ, nx, limbd, form='(I10,F8.4," - number of wavelengths, limb darkening")'
      for i=0L,nx-1 do printf, outQ, x(i), spq(i)
      printf, outQ, ' '+infile
      printf, outQ, fcblue*fcscale, fcred*fcscale, fc0*fcscale, $
              form= '(" FCblue=",E12.6," FCred=",E12.6," FC(0)=",E12.6)'
      free_lun, outQ
      printf, outU, nx, limbd, form='(I10,F8.4," - number of wavelengths, limb darkening")'
      for i=0L,nx-1 do printf, outU, x(i), spu(i)
      printf, outU, ' '+infile
      printf, outU, fcblue*fcscale, fcred*fcscale, fc0*fcscale, $
              form= '(" FCblue=",E12.6," FCred=",E12.6," FC(0)=",E12.6)'
      free_lun, outU
      printf, outV, nx, limbd, form='(I10,F8.4," - number of wavelengths, limb darkening")'
      for i=0L,nx-1 do printf, outV, x(i), spv(i)
      printf, outV, ' '+infile
      printf, outV, fcblue*fcscale, fcred*fcscale, fc0*fcscale, $
              form= '(" FCblue=",E12.6," FCred=",E12.6," FC(0)=",E12.6)'
      free_lun, outV
    endif
  endif
  
end