PRO shobwc1
  ;*** Display blue-whale calls.
  common control,ibc,ienv,nenviter,ifilt,iprint
  common param,nptfsig,nptenv,noffkern,envnaver,f1bh1,f2bh1
  common scandat,scanfile,lunscan,nbc,fnambc,iyearbc,daybc,dbbc,cenbc, $
    sigbc,thnsigbc,thfxbc
  common bsg,bfile,bsig,bsigfil,d0bf,ipksig,icensig,i1sig,i2sig
  common bnv,benv,ecen
  common dat2fit,npt,t,e,b
  str1=''
  ;*** Initialization (includes reading in scan list).
  init
  ;*** Loop over signals from scan list.
  FOR ibc=0,nbc-1 DO BEGIN
    print,ibc,dbbc(ibc)
    ;*** There is an interesting test event on r9125510.41b, at point
    ;*** number 1,303,100, which translates into day 255.45970858, or
    ;*** 11:01:58.82 on day 255.
    fnamtst=strtrim(fnambc(ibc))
    fnamtst=strmid(fnamtst,strlen(fnamtst)-12,12)
    daytst=daybc(ibc)/86400000.
    devent=255.4597086d; time the test event was found, using f0-kern = 16.0 s.
    dttest=.00011574d; ten seconds, in days
    test = ibc eq 619
    test = fnamtst eq 'r9125510.40c' and daytst gt (devent-dttest) and $
      daytst lt (devent + dttest); special event
    test = dbbc(ibc) ge 4.
    test=285.6458 lt daytst and 285.7708 gt daytst; fine series of calls
    IF test THEN BEGIN
;    IF ibc eq 619 THEN BEGIN
;    IF dbbc(ibc) gt 3. THEN BEGIN
;    IF fnamtst eq 'r9125510.40c' and daytst gt (devent-dttest) and $
;      daytst lt (devent + dttest) THEN BEGIN
      get1sig; Read in the time series, into bsig.
      filtsig; Filter the signal bsig, into bsigfil.
      IF ifilt gt 0 THEN bsig=bsigfil
      getenv; Calculate the envelope.
;      fitsig
;      read,'Press any key to make graph.',str1
      show1; plot something
;      IF iprint ge 1 THEN read,'Press any key to  continue.',str1
    ENDIF
  ENDFOR
  stop
END


PRO show1
  common control,ibc,ienv,nenviter,ifilt,iprint
  common param,nptfsig,nptenv,noffkern,envnaver,f1bh1,f2bh1
  common scandat,scanfile,lunscan,nbc,fnambc,iyearbc,daybc,dbbc,cenbc, $
    sigbc,thnsigbc,thfxbc
  common bsg,bfile,bsig,bsigfil,d0bf,ipksig,icensig,i1sig,i2sig
  common bnv,benv,ecen
  common dat2fit,npt,t,e,b
  str1=''; input buffer
  nptbsig=n_elements(bsig)
  ncen=long(ecen)
  tsig=lindgen(nptbsig)-ncen; time in ms, relative to center
  ngate=long(1000); Take this many points on either side of center.
  ngate=long(100); Take this many points on either side of center.
  i1=0>(ncen-ngate); first point to take
  i2=(ncen+ngate)<(nptbsig-1); last point to take
  t=tsig[i1:i2]
  b=bsig[i1:i2]
  e=benv[i1:i2]
  npt=n_elements(b)
  plot,t,b, $
    title=strmid(bfile,strlen(bfile)-12,12)+'  day '+ $
    string( long(daybc(ibc)/86400000.),format='(i3)' )+ ' '+ $
    string( long((daybc(ibc)/3600000.) mod 24.), format='(i2)' )+':' + $
    string( long((daybc(ibc)/60000.) mod 60.),format='(i2.2)')+':' + $
    string( (daybc(ibc)/1000.) mod 60., format='(f6.3)')+'  dbbc = '+ $
    string(dbbc(ibc),format='(f7.3)'), $
    xtitle='time (msec)', $
    thick=2
  ;*** Do a fit to first and third harmonics.
  ;*** a = [phi0, f0, phir, ah1, ah3]
  f0gs=16.; start the fit here
  phiguess,ngate,f0gs,phi0gs; estimate the phase
  phirgs=1.; arbitrary guess
  ah1=benv(ncen); Take amplitude guess from the envelope.
  ah3=ah1/2.

  a=[phi0gs,f0gs,phirgs,ah1,ah3]
  w=replicate(1.d,n_elements(b))
  bfit=curvefit(t,b,w,a,chisq=csq,tol=1.e-6,/noderivative, $
    function_name='funct2h')
  oplot,t,bfit,color='ff'x,thick=3
  phirel=a[2]
  IF a[3]*a[4] lt 0 THEN phirel=phirel+!dpi
  phirel=(phirel+20.*!dpi) mod (2.*!dpi); Force between 0 and 2 pi.
  xyouts,.15,.9,/normal,'f0:'+string(a[1],format='(f7.3)')+ $
    '  phi-rel:'+string(phirel,format='(f6.3)')+ $
    '  ah3/ah1:'+string(a[4]/a[3],format='(f9.4)')
  str1='y'
;  read,'Write jpeg file? ("y" or return) ',str1
  IF str1 eq 'y' THEN BEGIN
    jpegname=strmid(bfile,strlen(bfile)-10,4) + $
      string( long((daybc(ibc)/3600000.) mod 24.), format='(i2)' )+ $
      string( long((daybc(ibc)/60000.) mod 60.),format='(i2.2)')+ $
      string( long((daybc(ibc)/1000.) mod 60.), format='(i2.2)')+'.jpg'
    print,'jpegname: ',jpegname
    write_jpeg,jpegname,tvrd(true=1),true=1
  ENDIF
END


PRO funct2h, t, a, f
  ;*** Calculate blue-whale "B"-call function, harmonics 1 and 3.
  ;*** The theoretical function is 
  ;***    f[i]=a1*cos(phi1[i])+a3*cos(phi3[i]),
  ;*** where
  ;***    phi1[i]=phi0 + 2 pi f0 t[i]
  ;*** and
  ;***    phi3[i]=phi0 + phir + 2 pi 3*f0 t[i]
  ;*** a = [phi0,f0,phir,a1,a3]
  ;*** t is the array of times relative to the center of the call, in msec.
  ;***    phi[i] = phi0 + 2 pi f0 (t-t0) + 1/2 2 pi alpha (t-t0)^2 ,
  common control,ibc,ienv,nenviter,ifilt,iprint
  common dat2fit,npt,tpt,e,b
  common bnv,benv,ecen
  str1=''; input buffer
  phi0=a[0]; phase of first harmonic at the center of the waveform
  f0=a[1]; frequency of the first harmonic
  phir=a[2]; phase of second harmonic relative to first harmonic
  a1=a[3]; amplitude of first harmonic
  a3=a[4]; amplitude of third harmonic
  tsec=t/1000.
  phi1=phi0+2.*!dpi*f0*tsec; vector of phases, first harmonic 
  phi3=phir+3.*phi0+2.*!dpi*3.*f0*tsec; vector of phases, second harmonic
  f = a1*cos(phi1) + a3*cos(phi3)
  IF iprint ge 4 THEN print,a,total((f-b)^2)/(n_elements(t)-3), $
    format='(3f8.4,2f8.1,e12.4)'
  IF iprint ge 4 THEN oplot,t,f,color='ff0000'x
  IF iprint ge 6 THEN read,'Press any key to continue.',str1
END


PRO init
  common control,ibc,ienv,nenviter,ifilt,iprint
  common scandat,scanfile,lunscan,nbc,fnambc,iyearbc,daybc,dbbc,cenbc, $
    sigbc,thnsigbc,thfxbc
  common param,nptfsig,nptenv,noffkern,envnaver,f1bh1,f2bh1
  str1=''; input buffer
  device,decomposed=1; Select true color.
  plot,indgen(100)
  !p.background='ffffff'x; white background
  !p.color=0; default color is black
  plot,indgen(100)
  scanfile='/usr/data/bland03/pioneer/2091/244/r9124413.bwc.21217'
  scanfile='/usr/data/bland03/pioneer/2091/244/r9124413.bwc.21113'
  read,'scanfile '+scanfile+': ',str1
  IF strlen(str1) gt 0 THEN scanfile=str1

  ;*** ienv = 0 -> energy-weighted average
  ;*** ienv = 1 -> (energy)^(1.2)-weighted average
  ;*** ienv = 2 -> Gaussian fit
  ienv=1
  read,'ienv     '+string(ienv,format='(i6)')+': ',str1
  IF strlen(str1) gt 0 THEN ienv=fix(str1)

  ;*** nenviter = # of times to iterate selection of center of envelope
  nenviter=5

  ;*** ifilt = 0 -> no filtering of signal
  ;*** ifilt = 1 -> band-pass around first "B" harmonic
  ifilt=0
  read,'ifilt    '+string(ifilt,format='(i6)')+': ',str1
  IF strlen(str1) gt 0 THEN ifilt=fix(str1)

  ;*** iprint controls level of plotting and printing
  iprint=1
  read,'iprint   '+string(iprint,format='(i6)')+': ',str1
  IF strlen(str1) gt 0 THEN iprint=fix(str1)

  ;*** Number of time-series points in a full "B"-call signal.
  nptfsig=long(60000)
  nptenv=long(40000)
  noffkern=long(-5000)
  envnaver=500

  f1bh1=7.; lower frequency, B call, first harmonic
  read,'f1bh1    '+string(f1bh1,format='(f6.2)')+': ',str1
  IF strlen(str1) gt 0 THEN f1bh1=float(str1)

  f2bh1=20.; upper frequency, B call, first harmonic
  read,'f2bh1    '+string(f2bh1,format='(f6.2)')+': ',str1
  IF strlen(str1) gt 0 THEN f2bha=float(str1)
  getscdat
END
  

PRO getscdat
  ;*** Read in the information about calls found by the matched filter.
  common scandat,scanfile,lunscan,nbc,fnambc,iyearbc,daybc,dbbc,cenbc, $
    sigbc,thnsigbc,thfxbc
  openr,lunscan,scanfile,/get_lun,/swap_if_big_endian
  nbc=long(0)
  readu,lunscan,nbc
  print,'File ',scanfile,', ',strtrim(string(nbc),2),' events.'
  fnambc  = replicate('                                                  ',nbc)
  iyearbc = replicate(fix(0),   nbc)
  daybc   = replicate(double(0),nbc)
  dbbc    = replicate(float(0), nbc)
  cenbc   = replicate(float(0), nbc)
  sigbc   = replicate(float(0), nbc)
  thnsigbc= replicate(float(0), nbc)
  thfxbc  = replicate(float(0), nbc)

  readu,lunscan,fnambc,iyearbc,daybc,dbbc,cenbc,sigbc,thnsigbc,thfxbc
  close,lunscan
  free_lun,lunscan
END


PRO get1sig
  common control,ibc,ienv,nenviter,ifilt,iprint
  common param,nptfsig,nptenv,noffkern,envnaver,f1bh1,f2bh1
  common scandat,scanfile,lunscan,nbc,fnambc,iyearbc,daybc,dbbc,cenbc, $
    sigbc,thnsigbc,thfxbc
  common bsg,bfile,bsig,bsigfil,d0bf,ipksig,icensig,i1sig,i2sig
  str1=''
  ;*** Read in one blue-whale signal.
  bfile=strtrim(fnambc(ibc))
  strput,bfile,'b',strlen(bfile)-1
  openr,lunsig,bfile,/get_lun
  d0bf=double(0)
  npt=long(0)
  readu,lunsig,d0bf,npt
  fullsig=intarr(npt)
  readu,lunsig,fullsig
  close,lunsig
  free_lun,lunsig

  ;*** Select the segment containing the "B" call.
  ipksig=long(daybc(ibc)-d0bf); peak of matched-filter output
  icensig=ipksig-noffkern; subtract the kernel offset

  ;*** i1 and i2 are first and last points of the signal.
  i1sig=0>(icensig-nptfsig/2)<(npt-1)
  i2sig=(i1sig+nptfsig-1)<(npt-1)
  IF iprint ge 2 THEN print,'ipksig, icensig, i1sig, i2sig = ', $
    ipksig,icensig,i1sig,i2sig,format='(a,4i8)'
  bsig=fullsig[i1sig:i2sig]
END


PRO filtsig
  common control,ibc,ienv,nenviter,ifilt,iprint
  common param,nptfsig,nptenv,noffkern,envnaver,f1bh1,f2bh1
  common scandat,scanfile,lunscan,nbc,fnambc,iyearbc,daybc,dbbc,cenbc, $
    sigbc,thnsigbc,thfxbc
  common bsg,bfile,bsig,bsigfil,d0bf,ipksig,icensig,i1sig,i2sig
  str1=''
  zbsig=fft(bsig,-1)
  df=1000./nptfsig; frequency  step
  ;*** limits of frequency band
  i1=long(f1bh1/df)<(nptfsig-1)
  i2=long(f2bh1/df)<(nptfsig-1)
  IF f1bh1 gt 0. THEN zbsig(0)=0.; DC component
  IF i1 gt 1 THEN BEGIN
    zbsig(1:i1-1)=0.
    zbsig(nptfsig-i1+1:nptfsig-1)=0.
  ENDIF
  IF i2 eq nptfsig/2-1 THEN zbsig(nptfsig/2)=0.
  IF i2 lt nptfsig/2-1 THEN BEGIN
    zbsig(i2+1:nptfsig/2-1)=0.
    zbsig(nptfsig/2+1:nptfsig-i2-1)=0.
  ENDIF
  bsigfil=fix(fft(zbsig,1))
END


PRO getenv
  ;*** Calculate the envelope 
  common control,ibc,ienv,nenviter,ifilt,iprint
  common param,nptfsig,nptenv,noffkern,envnaver,f1bh1,f2bh1
  common scandat,scanfile,lunscan,nbc,fnambc,iyearbc,daybc,dbbc,cenbc, $
    sigbc,thnsigbc,thfxbc
  common bsg,bfile,bsig,bsigfil,d0bf,ipksig,icensig,i1sig,i2sig
  common bnv,benv,ecen
  str1=''

  IF iprint ge 3 THEN plot,bsig, $
    xtitle='one point per msec', $
    title=strmid(bfile,strlen(bfile)-12,12) + ': ipksig = ' + $
    strtrim(string(ipksig),2) + ', dbbc = ' + string(dbbc(ibc),format='(f6.3)')
  IF iprint ge 4 THEN read,'Press any key to continue',str1

  bsigsq=float(bsig)^2
  IF iprint ge 3 THEN plot,bsigsq, $
    xtitle='one point per msec', $
    title=strmid(bfile,strlen(bfile)-12,12) + ': ipksig = ' + $
    strtrim(string(ipksig),2) + ', dbbc = ' + string(dbbc(ibc),format='(f6.3)')
  
  ;*** Calculate an envelope by averaging envnaver points.
  envnaver=300
  ;*** For a sin-squared signal, the envelope is twice the average.
  esigsq=2.*smooth(bsigsq,envnaver,/edge_truncate)
  esig=sqrt(esigsq)
  IF iprint ge 3 THEN oplot,esigsq,color='00ff00'x
  IF iprint ge 4 THEN read,'Press any key to continue',str1

  ;*** Calculate the centroid of the envelope.  Weight with energy (== signal
  ;*** squared) for itype=0, and with sqrt(energy) for itype = 1.
  iptvect=lindgen(nptfsig); point numbers to go with esig
  pcen=fltarr(nenviter,2)
  FOR itype=0,1 DO BEGIN
    ;*** Start with the full signal length, then restrict to a central portion.
    i1=long(0)
    i2=nptfsig-1
    bs=esigsq
    IF itype eq 1 THEN bs = sqrt(esigsq)
    FOR ienviter=0,nenviter-1 DO BEGIN
      ;*** Calculate the energy-weighted centroid of the signal.
      pcen(ienviter,itype)=total(iptvect(i1:i2)*bs(i1:i2))/total(bs(i1:i2))
      ;*** Choose a signal symmetrical about the energy-weighted centroid, and
      ;*** realculate the energy-weighted centroid of the signal.
      ipcen=long(pcen(ienviter,itype))
      i1=ipcen-nptenv/2
      i2=i1+nptenv-1
     ENDFOR
  ENDFOR

  a=[1.,1.,1.,1.]; In case we skip the Gaussian fit.
  IF ienv eq 2 THEN BEGIN
    ;*** Carry out a Gaussian fit to the squared signal.
    t=lindgen(nptfsig); point numbers
    b=esigsq
    w=replicate(1.d,n_elements(b))
    backgdgs=total(b)/n_elements(b)/10.; initial guess
    signalgs=total(b)
    sigmags=5000.
    cengs=pcen(nenviter-1,0)
    a=[backgdgs,signalgs,cengs,sigmags]; vector of parameter initial guesses
    IF iprint ge 3 THEN plot,t,b,charsize=1.8,/nodata, $
      title='Gaussian fit to envelope',xtitle='point number'
    IF iprint ge 3 THEN oplot,t,b,color='00aaaa'x
    bfit=curvefit(t,b,w,a,chisq=csq,tol=1.e-6,/noderivative, $
      function_name='gauss1')
    IF iprint ge 3 THEN print,'a, csq = ',a,csq
    IF iprint ge 3 THEN oplot,t,bfit,thick=4,color=255
    IF iprint ge 3 THEN oplot,t,bfit-b,thick=2
    IF iprint ge 4 THEN read,'Press any key to continue.',str1
  ENDIF
  
  IF iprint ge 2 THEN print,'pcen = ',pcen
  IF  iprint ge 1 THEN print,'cen(E), cen(E^1/2), cen(G) = ', $
    pcen[nenviter-1,0],pcen[nenviter-1,1],a[2],format='(a,3f8.1)'

  ;*** Envelope and center to use in later fit.
  benv=esig
  ecen=pcen[nenviter-1,1<ienv]
  IF ienv eq 2 THEN ecen=a[2]

  dpeak=(d0bf+ipksig)/86400000.
  dpeak=daybc(ibc)/86400000.
  dcen=dpeak+(ecen-nptfsig/2-noffkern)/86400000.
  IF iprint ge 2 THEN plot,bsig, $
    xtitle='one point per msec', $
    title=strmid(bfile,strlen(bfile)-12,12) + ': dpeak=' + $
    string(dpeak,format='(f12.8)') + ', dcen=' + $
    string(dcen,format='(f12.8)') + $
    ', dbbc=' + string(dbbc(ibc),format='(f6.3)')
  IF iprint ge 2 THEN oplot,esig,color='00ff00'x
  IF iprint ge 2 THEN oplot,-esig,color='00ff00'x
  IF iprint ge 3 THEN read,'Press any key to continue',str1
END


PRO gauss1, t, par, f
  ;*** Fit to a Gaussian plus a constant background.
  ;*** The theoretical function is 
  ;***    f[i]= a + b*exp(-(t-t0)^2/(2 sig^2))/(sqrt(2 pi) sig)
  ;*** where a is the amplitude of the background, b is the amplitued of
  ;*** the signal, t0 is the center of the Gaussian, and sig is the width.
  common control,ibc,ienv,nenviter,ifilt,iprint
  a=par[0]
  b=par[1]
  t0=par(2)
  sig=par(3)
  f=a+b*exp(-(t-t0)^2/2./sig^2)/sqrt(2.*!pi)/sig
  IF iprint ge 4 THEN print,par
  IF iprint ge 4 THEN oplot,t,f,color='ff0000'x
END


PRO fitsig
  ;*** Fit the signal to a linear chirp.  A series of fits over increasing
  ;*** signal duration is used.
  common control,ibc,ienv,nenviter,ifilt,iprint
  common param,nptfsig,nptenv,noffkern,envnaver,f1bh1,f2bh1
  common scandat,scanfile,lunscan,nbc,fnambc,iyearbc,daybc,dbbc,cenbc, $
    sigbc,thnsigbc,thfxbc
  common bsg,bfile,bsig,bsigfil,d0bf,ipksig,icensig,i1sig,i2sig
  common bnv,benv,ecen
  str1=''; input buffer

  ;*** Display the signal
  ntot=nptfsig
  tsig=findgen(ntot)/1000.; times in seconds
  IF iprint ge 1 THEN plot,tsig,bsig,charsize=1.8,/nodata, $
    title='Blue-Whale "B" Call', $
    subtitle=strmid(bfile,strlen(bfile)-12,12)+ $
    '  '+string(long((d0bf+ipksig)/3600000.) mod 24,format='(i2)')+ $
    ':'+string(long((d0bf+ipksig)/60000.) mod 60,format='(i2.2)')+ $
    ':'+string(((d0bf+ipksig)/1000.) mod 60.,format='(f6.3)'), $
    xtitle='time (sec)', $
    ytitle='signal (ADC units)'
  IF iprint ge 1 THEN oplot,tsig,bsig,color='00aaaa'x
  IF iprint ge 2 THEN oplot,tsig,benv,thick=4,color=255
  IF iprint ge 2 THEN oplot,tsig,-benv,thick=4,color=255
  IF iprint ge 2 THEN read,'Press any key to continue.',str1

  ;*** a = [phi0, f0, dfdt]
  f0gs=16.
  ngate=long(100)
  phiguess,ngate,f0gs,phi0gs
  a=[phi0gs,16.,-.05]
  linchrp,100,a
  linchrp,300,a
  linchrp,500,a
  linchrp,1000,a
  linchrp,3000,a
END


PRO phiguess,ngate,f0gs,phi0gs
  ;*** Estimate the signal phase, for the initial guess in the lsq fit.
  common bsg,bfile,bsig,bsigfil,d0bf,ipksig,icensig,i1sig,i2sig
  common bnv,benv,ecen
  str1=''; input buffer
  ncen=long(ecen)
  tsig=lindgen(n_elements(bsig))-ncen; t-t0, in point numbers
  i1=0>(ncen-ngate); first point
  i2=(ncen+ngate)<(n_elements(bsig)-1); last point
  t=tsig[i1:i2]
  b=bsig[i1:i2]
  npt=n_elements(b)
  omega=f0gs/1000.*2.*!pi; omega in radians per ms
  c=total(b*cos(omega*t))
  s=total(b*sin(omega*t))
  c1=total(cos(omega*t)^2)
  c2=total(sin(omega*t)*cos(omega*t))
  s1=total(cos(omega*t)*sin(omega*t))
  s2=total(sin(omega*t)^2)

  d=s2*c1-c2*s1
  cosphi0=(s2*c-c2*s)/d
  sinphi0=(s1*c-c1*s)/d
  phi0gs=atan(sinphi0,cosphi0)
END


PRO linchrp,ngate,a
  ;*** Carry out a linear-chirp least-squares fit to the signal.
  common control,ibc,ienv,nenviter,ifilt,iprint
  common param,nptfsig,nptenv,noffkern,envnaver,f1bh1,f2bh1
  common scandat,scanfile,lunscan,nbc,fnambc,iyearbc,daybc,dbbc,cenbc, $
    sigbc,thnsigbc,thfxbc
  common bsg,bfile,bsig,bsigfil,d0bf,ipksig,icensig,i1sig,i2sig
  common bnv,benv,ecen
  common dat2fit,npt,t,e,b
  str1=''; input buffer
  ntot=nptfsig; number of signal points
  tsig=findgen(ntot)/1000.; times in seconds

  npt=2*ngate+1
  t=tsig(0>(ecen-ngate):(ecen+ngate)<(ntot-1))
  b=bsig(0>(ecen-ngate):(ecen+ngate)<(ntot-1))
  e=benv(0>(ecen-ngate):(ecen+ngate)<(ntot-1))

  IF iprint ge 3 THEN plot,t,b,charsize=1.8,/nodata, $
    title='Selected Region of Whale Call',xtitle='time (sec)', $
    ytitle='signal (ADC units)'
  IF iprint ge 3 THEN oplot,t,b,color='00aaaa'x
  IF iprint ge 3 THEN oplot,t,e,thick=4,color=255
  IF iprint ge 3 THEN oplot,t,-e,thick=4,color=255
  IF iprint ge 4 THEN read,'Press any key to continue.',str1

  w=replicate(1.d,n_elements(b))
  bfit=curvefit(t,b,w,a,chisq=csq,tol=1.e-6,/noderivative, $
    function_name='functlch')
  IF iprint ge 3 THEN oplot,t,bfit-b,thick=2
  IF iprint ge 1 THEN print,ngate,a,csq,format='(i6,f9.4,2f10.4,f12.2)'
  IF iprint ge 4 THEN read,'Press any key to continue.',str1
END


PRO functlch, t, a, f
  ;*** Calculate blue-whale "B"-call function for fitting program.
  ;*** The theoretical function is 
  ;***    f[i]=e[i]*cos(phi[i],
  ;*** where e is the envelope of the signal being fitted, and the phase
  ;*** phi is given by
  ;***    phi[i] = phi0 + 2 pi f0 (t-t0) + 1/2 2 pi alpha (t-t0)^2 ,
  ;*** where t0 is the time of the center of the whale call, as estimated
  ;*** from the envelope.  The envelope and the index of its center are
  ;***  stored in the common block bwcdata.
  common control,ibc,ienv,nenviter,ifilt,iprint
  common dat2fit,npt,tpt,e,b
  common bnv,benv,ecen
;  common param,nptfsig,nptenv,noffkern,envnaver,f1bh1,f2bh1
;  common scandat,scanfile,lunscan,nbc,fnambc,iyearbc,daybc,dbbc,cenbc, $
;    sigbc,thnsigbc,thfxbc
;  common bsg,bfile,bsig,bsigfil,d0bf,ipksig,icensig,i1sig,i2sig
  str1=''; input buffer
  t0=ecen/1000.d; one point per millisecond
  phi0=a[0]; phase at t0
  f0=a[1]; frequency at t0
  alpha=a[2]; frequency sweep rate
  phi=phi0+2.*!dpi*f0*(t-t0)+0.5*2.*!dpi*alpha*(t-t0)^2
  f=e*cos(phi)
  IF iprint ge 4 THEN print,a,total((f-b)^2)/(n_elements(t)-3), $
    format='(3f12.8,e12.4)'
  IF iprint ge 4 THEN oplot,t,f,color='ff0000'x
  IF iprint ge 6 THEN read,'Press any key to continue.',str1
END