PRO plotbw2
  ;*** Plots for baby-B whale, Feb. 26, 2005, RWB.
  ;*** Make some plots for Blue-Whale calls.  The data should be SAVEd by
  ;*** rdfit.pro .  The variables which should be in memory are
  ;***	npt:	number of points in vectors
  ;***	f0	central frequency
  ;***	alpha	negative of slope
  ;***	phi0	phase at center of call
  ;***	res2	residuals per degree of freedom
  ;***	nchi2	res2^2/(sig*fit)
  ;***	index	number of event in input scan list
  ;***  infile	input file, complete with path
  ;***  fname	file name, w/o path

  str1=''; input buffer
  plot,indgen(100)
  !p.color=0
  !p.background='ffffff'x
  restore

  ;******************************
  ;******* f0 versus alpha ******
  ;******************************
  makplot,f0,12.,.1,80,alpha,-40.,1.,80,fname+'  All Events', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2a.jpg'

  ;*** Same distributions, for a restricted region coresponding to the "outer
  ;*** hexagon."
  makplot,f0,15.5,.01,100,alpha,-.2,.005,100,fname+'  All, blowup', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2b.jpg'

  ;*** Same distributions, for a restricted region coresponding to the "inner
  ;*** hexagon."
  makplot,f0,15.7,.010,60,alpha,-.10,.005,60,fname+'  All, central', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2c.jpg'

  ;*** Same distributions, for the "inner sphere" region.
  makplot,f0,15.80,.005,80,alpha,0.0,.002,50,fname+'  All, in. sph', $
    'central frequency (Hz)','-dfdt (Hz/sec)',IGAUSS=1
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2d.jpg'


  ;************************************
  ;******* nchi2 versus res2 ********
  ;************************************
  ;*** all events
  makplot,nchi2,0.,.04,50,res2,0.,1.e8,100,fname+'  All Events', $
    'nchi2 (res^2/(|sig|*|fit|))','(squared residuals)/dof'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2e.jpg'

  ;*** blowup
  makplot,nchi2,0.,.02,50,res2,0.,1.e7,50,fname+'  All, blowup', $
    'nchi2 (res^2/(|sig|*|fit|))','(squared residuals)/dof'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2f.jpg'

  ;***********************************
  ;*******  Select events.  **********
  ;***********************************

  ;*** Convergent events - first five fits converged (except maybe the first
  ;*** fit), OR the cheating fit improved the chi-squared by a certain margin.
  itest = exits ge 30
  ix=where(itest)

  makplot,f0[ix],12.,.1,80,alpha[ix],-40.,1.,80,fname+'  Convergent', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2g.jpg'

  ;*** Same distributions, for a restricted region coresponding to the "outer
  ;*** hexagon."
  makplot,f0[ix],15.5,.01,100,alpha[ix],-.2,.005,100,fname+'  Cvg, blowup', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2h.jpg'

  ;*** Same distributions, for a central region coresponding to the "inner
  ;*** hexagon."
  makplot,f0[ix],15.7,.010,60,alpha[ix],-.10,.005,60,fname+'  Cvg, central', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2i.jpg'

  ;*** Good convergent events - convergent as above, and nchi2 < .7.
  itest = exits ge 30 and nchi2 lt .7
  ix=where(itest)

  ;*** Full range plot.
  makplot,f0[ix],12.,.1,80,alpha[ix],-40.,1.,80,fname+'  Good, cvgt', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2j.jpg'

  ;*** Blow up in alpha only.
  makplot,f0[ix],12.,.1,80,alpha[ix],-.5,.02,50,fname+'  Good, cvgt', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2k.jpg'

  ;*** Same distributions, for a restricted region coresponding to the "outer
  ;*** hexagon."
  makplot,f0[ix],15.5,.01,100,alpha[ix],-.2,.005,100,fname+'  GC, blowup', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2l.jpg'

  ;*** Same distributions, for a central region.
  makplot,f0[ix],15.7,.010,60,alpha[ix],-.10,.005,60,fname+'  GC, central', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2m.jpg'


  ;*** Start a series of alpha vs f0, central region, for different cuts
  ;*** to select events with accurately measured parameters.

  ;*** All convergent events.
  itest = exits ge 30
  ix=where(itest)
  makplot,f0[ix],15.7,.010,60,alpha[ix],-.10,.005,60,fname+'  Cvg', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2n.jpg'

  ;*** All convergent events, nchi2 gt 1..
  itest = exits ge 30 and nchi2 gt 1.
  ix=where(itest)
  makplot,f0[ix],15.7,.010,60,alpha[ix],-.10,.005,60,fname+'  Cvg, nc2>1.', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2o.jpg'

  ;*** All convergent events, nchi2 lt 1..
  itest = exits ge 30 and nchi2 lt 1.
  ix=where(itest)
  makplot,f0[ix],15.7,.010,60,alpha[ix],-.10,.005,60,fname+'  Cvg, nc2<1.', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2p.jpg'

  ;*** All convergent events, nchi2 lt .7.
  itest = exits ge 30 and nchi2 lt .7
  ix=where(itest)
  makplot,f0[ix],15.7,.010,60,alpha[ix],-.10,.005,60,fname+'  Cvg, nc2<.7', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2q.jpg'

  ;*** All convergent events, nchi2 lt .5.
  itest = exits ge 30 and nchi2 lt .5
  ix=where(itest)
  makplot,f0[ix],15.7,.010,60,alpha[ix],-.10,.005,60,fname+'  Cvg, nc2 <.5', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2r.jpg'

  ;*** All convergent events, nchi2 lt .3.
  itest = exits ge 30 and nchi2 lt .3
  ix=where(itest)
  makplot,f0[ix],15.7,.010,60,alpha[ix],-.10,.005,60,fname+'  Cvg, nc2 <.3', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2s.jpg'

  ;*** All convergent events, nchi2 lt .1.
  itest = exits ge 30 and nchi2 lt .1
  ix=where(itest)
  makplot,f0[ix],15.7,.010,60,alpha[ix],-.10,.005,60,fname+'  Cvg, nc2 <.1', $
    'central frequency (Hz)','-dfdt (Hz/sec)'
  read,'Write jpeg file?  ("y" ou merde) ',str1
  IF str1 eq 'y' THEN w1jpeg,'plotbw2t.jpg'

  dtest=cday(index)/86400000.+(iyear-1)*365.
  itest=exits ge 30 and dtest gt 385.5 and dtest lt 385.8

END


PRO makplot,x,xmin,xbinsize,xnbins,y,ymin,ybinsize,ynbins,title, $
  xtitle,ytitle,igauss=igauss
  window,xsize=800,ysize=600
  str1=''; input buffer
  !p.multi=[0,2,2,0,0]
  chsave=!p.charsize; Save current character size.
  !p.charsize=1.6; Set larger characters.

  xmax=xmin+xbinsize*xnbins
  ymax=ymin+ybinsize*ynbins
  plot,x,y,psym=3,xrange=[xmin,xmax],yrange=[ymin,ymax], $
    title=title,xtitle=xtitle,ytitle=ytitle,xstyle=1,ystyle=1
  ntot=n_elements(x)
  test1= x ge xmin and x lt xmax and y ge ymin and y lt ymax
  npt=long(total(test1))
  xyouts,/normal,.3,.93,string(ntot),charsize=1
  xyouts,/normal,.3,.91,string(npt), charsize=1

  ;*** If a gaussian fit is specified, carry it out.
  IF n_elements(igauss) ne 0 THEN BEGIN
    ibins=histogram(test1,reverse_indices=ix)
    xplot=x(ix(ix(1):ix(2)-1))
    yplot=y(ix(ix(1):ix(2)-1))
    gfit2d,xplot,yplot,rho,IPLOT=1
  ENDIF

  makhist,y[where(test1)],ymin,ybinsize,ynbins,'y projection',ytitle
  makhist,x[where(test1)],xmin,xbinsize,xnbins,'x projection',xtitle

  !p.multi=0; Restore single-plot mode.
  !p.charsize=chsave; Restore previous character size
END


PRO makhist,f,fmin,fbinsize,fnbins,ftitle,fxtitle
  str1=''; input buffer
  yhist=histogram(f,min=fmin,binsize=fbinsize,nbins=fnbins,reverse_indices=ix)
  nhist=n_elements(yhist)
  xval=fmin+findgen(fnbins)*fbinsize
  plot,xval,yhist,psym=10, title=ftitle,xtitle=fxtitle,xstyle=1,ystyle=1
  oplot,xval,yhist,psym=10,thick=1.5
  ;*** Calculate mean and standard deviation.
  fsel=f[ix[ix[0]:ix[fnbins]-1]]
  ntot=n_elements(f)
  nsel=ix[fnbins]-ix[0]
  fmean=mean(fsel)
  fstddev=stddev(fsel)

  ;*** Calculate a 1.5-sigma mean and stdev, where only events within 1.5 
  ;*** sigmas are used to recursively recalculate the mean and stdev
  ;*** (correcting the stdev for the missing tails, of course).
  npass=9
  fnout=1.5
  ;*** Start with values above.
  fmean1=fmean
  fstddev1=fstddev
;  print,'    fmean1  fstddev1'
;  print,fmean1,fstddev1
  FOR i=0,npass DO BEGIN
    ;*** Iterate, taking fnout sigmas out from the mean.
    test2=f gt fmean1-fstddev1*fnout and f lt fmean1+fstddev1*fnout
    fmean1=mean(f[where(test2)])
    fstddev1=stddev(f[where(test2)])
    ;*** Correct for only going out fnout sigmas in the calculation.
    falf=1.-2.*gauss_pdf(-fnout); prob. inside +/- fnout sigmas
    cfact=sqrt(  falf/(falf-2./sqrt(2.*!pi)*fnout*exp(-fnout^2/2.))  )
    fstddev1=fstddev1*cfact
;    print,fmean1,fstddev1
  ENDFOR
  oplot,xval,nsel/sqrt(2.*!pi)/fstddev1*fbinsize* $
    exp(-(fmin+findgen(nhist)*fbinsize-fmean1)^2/(2.*fstddev1^2))
;  read,'Press any key to continue.',str1
  

  xval=.35
  yval=.43
  IF !p.multi[0] eq 2 THEN xval=xval+.5
  IF !p.multi[0] eq 2 THEN yval=yval+.5
  xyouts,/normal,xval,yval,string(ntot),charsize=1
  xyouts,/normal,xval,yval-.02,string(nsel),charsize=1
  xyouts,/normal,xval,yval-.04,string(fmean),charsize=1
  xyouts,/normal,xval,yval-.06,string(fstddev),charsize=1
  xyouts,/normal,xval,yval-.08,string(fmean1),charsize=1
  xyouts,/normal,xval,yval-.10,string(fstddev1),charsize=1
END


PRO w1jpeg,fname
  str1=''; input buffer
  read,'fname '+fname+' :',str1
  IF strlen(str1) gt 0 THEN fname=str1
  write_jpeg,fname,tvrd()
END


PRO gfit2d,x,y,rho,IPLOT=iplot
  ;*** Calculate a 2-dimensional Gaussian fit to a sample of coordinates
  ;*** {(xi,yi)}.
  npt=n_elements(x)
  IF npt lt 3 THEN RETURN
  ;*** Everything becomes floating point from here on.
  x2 = total(x*x)/npt
  y2 = total(y*y)/npt
  xy = total(x*y)/npt
  rho = 1./(x2*y2-4.*xy^2)
  alpha = sqrt(rho)*y2
  beta = sqrt(rho)*x2
  gamma = sqrt(alpha*beta-1.)
END