	subroutine flake
c	makes file for zeta plotter using the bdms system
c	assumes de, elow and noe are same for all triggers
c	other than lwbias trigger which may have different values
c
	common/io/b(100,2),ee(100,2),ben(100,2),bl(20,2),be(20,2),
&	blen(20,2),n,ln,nfsq,inpar,thii,tloo,avtt,lwprob,hiunit
	logical hiunit,nfsq,lwprob,doglab,doxlab,doylab,dodots
c
	include '/u1/gcf/e350/crystal/cb/wenbin'
	include '/u1/gcf/e350/crystal/cb/foxin'
	include '/u1/gcf/e350/crystal/cb/clouds'
	include '/u1/rgk/icicle/cb/ztpltc'
	include '/u1/gcf/e350/crystal/cb/hotsun'
	include '/u1/gcf/e350/crystal/cb/where'
	include '/u1/gcf/e350/crystal/cb/it'
	include '/u1/gcf/e350/crystal/cb/igr'
	include '/u1/gcf/e350/crystal/cb/const'
	include '/u1/gcf/e350/crystal/cb/datbin'
	include '/u1/gcf/e350/crystal/cb/they'
	include '/u1/gcf/e350/crystal/cb/beam'
	include '/u1/gcf/e350/crystal/cb/parms'
	include '/u1/gcf/e350/crystal/cb/abssiz'
	include '/u1/gcf/e350/crystal/cb/ooflak'
	common/lbias/lbias(7)
	dimension bnew(180,2),enew(180,2),dotted(180,2),iebeg(7),ie1min(2)
	dimension nam1(6),nam2(2),thmax(4),iinia(2),acpi(150)
	dimension je2(7),je1(7),ilow(7),acth(4)
	equivalence (iinia,iini1),(iinia(2),iini2)
	logical dotted,start,a90,only1
	character*10 nam1,nam2,nono
c
	data nam1/7h$3p$u0-,7h$3p$u0+,7h$2k$u0-,7h$2k$u0+,4h$~1p,3h$1p/
	data nam2/7h$3p$u70,3h$3e/,nplot/0/,ccxl/265./,ccyl/100./
	data lwprob/.false./,labely/500/,uplimt/1.e38/,doglab/.true./
	data doylab/.true./,doxlab/.true./,dodots/.false./
c	dodots=.true.: use dotted lines for data points not in fit
c	not implemented at present.  extra dpt made squares appear - see p88
c
300	format('options=x;lim='e10.3,',',e10.3,';ticks='i10,','i10';'/
&	'options.=y;lim='e10.3','e10.3';loglin=log 'i3';')
301	format('pltype=xlabel;'/
&	'pltext=.       E (GeV)         .;')
c302	format(50hylabel    $1$d$u72$3s$1/dtdx mb/$2(g$1e$2v/$1c)$u72)
302	format('pltype=ylabel;'/
&	'pltext=. cross section.;')
c502	format('ylabel     $1d$u02$3s$1/dtd$2e $3m$1b/$2g$1e$2v$u03'/
c&	'labelspecsplot$lc',7x,'4',4x,'1 -75.',5x,f10.3,f10.3,'200.')
c303	format(10hglabel    ,a7,6h$1p$5(,2a7)
c305	format(12hsymbol    1.)
306	format('dpt=')
307	format('*e')
309	format(2x,e15.5,'+-',e15.5,','e15.5,'+-',e15.5,';')
310	format(2x,e15.5,',',e15.5,';')
311	format('style='i2';'/'pt=')
c313	format(4hccxy,6x,4f10.3)
c314	format(6hdlabel,8x,1h4,4x,1h0,f10.3,7x,i3,2h$-,i1,2x,3a10)
315	format('family.=plot'i2';'/'plot=hippo;')
316	format('add')
c317	format(20hsymbol     102000001)
c321	format(25hglblspecs             4.0)
c323	format(6hdlabel,8x,i1,4x,1h0,2f10.3,f4.2,8h<-$1t$2<,f4.2,1h ,
c&	'$2g$1e$2v$u02')
323	format('pltype=title;'/
&	'pltext='f4.2'<-t<'f4.2' GeV**2;')
325	format('pltype=title;'/
&	'pltext=.'f3.2'<-t<'f3.2' GeV**2.;')
326	format('pltype=glabel;'/
&	'pltext=.'f3.2'<-t<'f4.2' GeV**2.;')
c324	format(6hdlabel,8x,i1,4x,1h0,2f10.3,6h$3a$2=,f6.3,
c&	5h$5p$2,f5.3)
c327	format(6hdlabel,8x,i1,4x,1h0,2f10.3,7h$3a$2=0,f4.3,6h$5p$20,
c&	f4.3)
c328	format(6hdlabel,8x,i1,4x,1h0,2f10.3,8h$3a$2=-0,f4.3,6h$5p$20,
c&	f4.3)
330	format(' maximum y value = ',e12.4,' top of plot = ',e12.4)
c331	format(6harrows/11hctype     0/5hcurve/10x,2e15.5/10x,2e15.5/
c&	8hendcurve/5hcurve,2(/10x,2e15.5)/8hendcurve/8hnoarrows)
c334	format(8haxislbls,10x,2hon,8x,2hon)
c339	format(6hdlabel,8x,i1,4x,1h0,2f10.3,1h0,f3.2,7h<$1t$2<,1h0,f3.2,
c&	1h ,19h$2(g$1e$2v/$1c)$u72)
c340	format(6hdlabel,8x,1h4,4x,1h0,f10.3,7x,i3,2h$-,i3,2x,3a10)
c341	format(8hboundary,8x,4haxis,6x,4hline,6x,4haxis,6x,4haxis)
c342	format(9hylblspecs,11x,3h4.0)
c343	format('glabel')
c344	format('xlabel'/8haxislbls,9x,3hoff,7x,2hon)
c346	format(9hylblspecs,11x,3h4.0,5h 63.0)
c347	format(6hdlabel,8x,i1,4x,1h0,2f10.3,6h$3a$2=,f5.2,6h$5p$20,f3.2)
c348	format(6hdlabel,8x,i1,4x,1h0,2f10.3,7h$3a$2=0,f3.2,6h$5p$20,f3.2)
c349	format(6hdlabel,8x,i1,4x,1h0,2f10.3,8h$3a$2=-0,f3.2,6h$5p$20,f3.2)
350	format('dptext=1;'/'dpparm=0.01;'/'dpt=')
c
403	format(' no plot produced because set',i2,' has bin size',f6.2,
&	' different from desired',f6.2,' lwbias is set no',i2,
&	' group no',i2,' igrp=',7i2)
404	format(' cross sections are ',6f8.4)
c
	if((.not.doplot).and.(.not.dolist)) return
	if (nplot.eq.0) goto 93
94	nplot=nplot+1
c	write(3,334)
	do 86 j=1,2
	iinia(j)=0
	do 87 l=1,ntrdis
	if (igrp(l).eq.j.and.(.not.empty(l))) goto 88
87	continue
	goto 86
88	iinia(j)=l
86	continue
	if(iini1.eq.0.and.iini2.eq.0)return
	if(.not. dolist) goto 55
	if(iinia(2).ne. 0) goto 55
	if (nfs) goto 55
	if (de(iinia(1)).ne.de(iinia(2))) call abort
	if((lwbias(1).eq.0).or.(lwbias(2).eq.0)) goto 55
	if(de(lwbias(1)).ne. de(lwbias(2))) call abort
55	do 16 l=1,ntrdis
	if (l.eq.lwbias(1).or.l.eq.lwbias(2)) goto 85
	iini=iinia(igrp(l))
	iebeg(l)=max0(1,ifix((ebeg(l)-elow(iini))/de(iini)-0.000001)+2)
	je1(l)=ifix((e1(l)-elow(iini))/de(iini)+1)
	je2(l)=ifix((e2(l)-elow(iini))/de(iini)+1)
	ilow(l)=(elow(l)-elow(iini))/de(iini)
	goto 16
85	iebeg(l)=max0(1,ifix((ebeg(l)-elow(l))/de(l)-0.000001)+2)
	je1(l)=ie1(l)
	je2(l)=ie2(l)
16	continue
	t1=abs(tlo)
	t2=abs(thi)
	lgplot=3
	if (t1.ge.0.4) lgplot=4
	ichi=1
	if (t1.ge.2.0) ichi=2
	eplx=100.
	ephx=210.
	ego=enhigh((nobeam+1)/2)
c	use energy as variable
	ego=1.0
	fudm=fdmicr*ego/abs(tlo-thi)
	ccyh=100.+183.3*float(lgplot)
	ccxlus=ccxl
	ccxhus=815.
	if (iglow) ccxlus=448.
c	write(3,313) ccxlus,ccxhus,ccyl,ccyh
c	write(3,321)
	write(3,315) nplot
	if (doxlab) write(3,301)
c	if (.not.doxlab) write(3,344)
	hiunit=.false.
	if (nfs.or.t1.gt.0.99.or.(t1.ge.1.0.and.ipart.eq.2).or.ego.eq.1.0)
&	goto 29
	fudm=.001*fudm
	hiunit=.true.
c	if (ego.eq.1.0.or.dolist) write(3,502)
29	blm=0.
c	average over data sets and set up y limits of plot, making limits same
c	for the 2 groups
	do 92 j=1,2
	iini=iinia(j)
	if (iini.eq.0) goto 92
c	check that binning is same
	dd=de(iini)
	do 39 l=1,ntrdis
	if (igrp(l).ne.j) goto 39
	if (l.eq.lwbias(j)) goto 39
	if (de(l).eq.dd) goto 39
	print 403,l,de(l),dd,lwbias(j),j,(igrp(lj),lj=1,ntrdis)
	iinia(j)=0
	goto 92
39	continue
	noet=noe(iini)
	istrt=max0(1,ifix((eplx-elow(iini))/dd)+1)
	if (noet-istrt+1.gt.180) call abort
	do 2 i=istrt,noet
	dotted(i,j)=.true.
	dd=0.
	bb=0.
	do 9 l=1,ntrdis
	if (igrp(l).ne.j) goto 9
	if (l.eq.lwbias(j)) goto 9
	if (i.lt.iebeg(l)) goto 9
	if ((.not.dodots).or.(dodots.and.i.le.je2(l).and.i.ge.je1(l)))
&	dotted(i,j)=.false.
	lj=i-ilow(l)
	cor=1.
	if (nfs) cor=core(l)
	eb2=(de(l)*ebinb(lj,l)/(fudm*cor))**2
	bb=bb+binb(lj,l)*eb2*fudm*cor/de(l)
	dd=dd+eb2
9	continue
	if (dd.eq.0.) goto 10
	bnew(i,j)=bb/dd
	enew(i,j)=1./sqrt(dd)
	blm=amax1(blm,bnew(i,j)+enew(i,j))
	goto 24
10	bnew(i,j)=1.e-20
	enew(i,j)=1.e-20
24	continue
2	continue
92	continue
c
	if (iini2.eq.0.and.iini1.eq.0) return
	do 25 j=1,2
	if (iinia(j).eq.0) goto 25
	ill=lwbias(j)
	if (ill.eq.0) goto 25
	e1min=1.e10
	do 11 l=1,ntrdis
	if (igrp(l).ne.j) goto 11
	if (l.eq.ill) goto 11
	e1min=amin1(e1min,ebeg(l))
11	continue
	del=de(ill)
	elowl=elow(ill)
	ie1min(j)=(e1min-elowl)/del
	ixx=iebeg(ill)
	if (ixx.gt.ie1min(j)) goto 25
	ixz=ie1min(j)
	if (iglow) goto 25
	do 21 i=ixx,ixz
	if (elowl+float(i)*del.lt.eplx) goto 21
	ebb=ebinb(i,ill)
	if (ebb.gt.0.) blm=amax1(blm,(binb(i,ill)+1./ebb)*fudm/del)
21	continue
25	continue
	if(dolist .and. .not. doplot) goto 42
	if(dolist) blm=blm/ego
	ipa=alog10(blm)+1
	epa=10.**ipa
	if (ipa.lt.0.) epa=epa*10.
	if (blm.gt.0.3*epa) goto 32
	epa=0.3*epa
32	print 330,blm,epa
	epa=amin1(uplimt,epa)
14	eqa=0.
	eqa=epa/(10.**lgplot)
	t4=t1+0.0000001
	if (t1.ge.1.0.and.t2.ge.1.0) write(3,323) t1,t2
	if (t1.lt.1.0.and.t2.ge.1.0) write(3,326) t1,t2
	if (t1.lt.1.0.and.t2.lt.1.0) write(3,325) t4,t2
	write(3,300) eplx,ephx,nlab,ntkx,eqa,epa,lgplot
c	if (eqa.lt.0.0099) write(3,346)
	plq=ccyh-ccyl+100.
	plr=-50.
	if (eqa.lt.0.0009) plr=-20.
c	if (doylab) write(3,502) plr,plq
c
c	initialize lwbias array
	do 56 lj=1,20
56	bl(lj,2)=-1
c
c	write out data
42	do 26 iq=1,2
	if (iinia(iq).eq.0) goto 26
	iini=iinia(iq)
	ddi=de(iini)
	istrt=max0(1,ifix((eplx-elow(iini))/ddi)+1)
	e=elow(iini)+(float(istrt)-0.5)*ddi
	noet=noe(iini)
	if(.not. doplot) goto 52
c	if (iq.eq.1) goto 54
c	ccxlus=25.0
c	ccxhus=575.0
c	if (lwbias(2).eq.0.or.iglow) ccxhus=392.0
c	write(3,313) ccxlus,ccxhus,ccyl,ccyh
c	if (lwbias(2).ne.0) goto 54
c	eplai=70.
c	write(3,300) eplai,ephx,eqa,epa
c	write out glabel header
54	if (.not.doglab) goto 67
c	if (ipart.eq.1) write(3,303) nam1(jpart(iq)),nam2(ipart),nono
	goto 52
67	continue
c	write(3,343)
52	nloop=iocon1(iq)
	ii=0
	start=.true.
	n=0
	if((iq.eq.1).or.(lwbias(2).ne.0)) ln=0
	a90=.false.
	do 3 i=istrt,noet
	if (enew(i,iq).lt.1.e-19) goto 3
	ii=ii+1
	if (ii.eq.1) goto 61
	if(iloop.ne.nloop) goto 62
61	eav=0.
	bbav=0.
	ebav=0.
	iloop=0
	range=0.
cveto	if (e.lt.180..or.a90) goto 62
cveto	nloop=ichi
cveto	a90=.true.
cveto	if (nfs.and.t1.lt.0.22) write(3,350)
c
62	iloop=iloop+1
	range=range+0.5*ddi
	eav=eav+e
	bbav=bbav+bnew(i,iq)
	ebav=ebav+enew(i,iq)**2
	if(iloop.ne.nloop.and.i.ne.noet) goto 3
c
	eav=eav/iloop
	bbav=amax1(eqa,bbav/iloop)
	ebav=sqrt(ebav)/iloop
	if (start) goto 18
	if (dotted(i,iq).and.(.not.dotted(i-nloop,iq))) goto 15
	if ((.not.dotted(i,iq)).and.dotted(i-nloop,iq)) goto 17
	goto 19
c15	write(3,317)
15	write(3,306)
	goto 19
c	write(3,305)
17	write(3,306)
	goto 19
c18	if (dotted(i,iq)) write(3,317)
c	if (.not.dotted(i,iq)) write(3,305)
18	write(3,306)
	start=.false.
19	if(dolist) goto 80
	range=range/ego
	eav=eav/ego
	write(3,309) eav,range,bbav,ebav
	goto 3
80	y=bbav/ego
	ye=ebav/ego
c   **iq=2 means pi-plus**
	n=n+1
	b(n,iq)=bbav
	ee(n,iq)=ebav
	ben(n,1)=eav-range
	ben(n,2)=eav+range
	write(3,309) eav,range,y,ye
3	e=e+ddi
c
c	write out lwbias data
	ill=lwbias(iq)
	if (ill.eq.0) goto 1
c	if (.not.dodots) write(3,305)
c	if (dodots) write(3,317)
	del=de(ill)
	ixx=iebeg(ill)
	e=elow(ill)+(ixx-.5)*del
	nloop=iocon2(iq)
	ixz=min0(ie1min(iq),noe(ill))
	if (ixx.gt.ixz) goto 1
	write(3,306)
	do 12 i=ixx,ixz
	if (i.eq.ixx) goto 71
	if(iloop.ne.nloop) goto 72
71	eav=0.
	bbav=0.
	ebav=0.
	iloop=0
	range=0.
c
72	iloop=iloop+1
	range=range+.5*del
	eav=eav+e
	cor=1.
	if (nfs) cor=core(ill)
	bbav=bbav+fudm*binb(i,ill)*cor/del
	if (ebinb(i,ill).gt.0.) ebav=ebav+(fudm*cor/(del*ebinb(i,ill)))**2
	if (iloop.ne.nloop.and.i.ne.ixz) goto 12
c
	eav=eav/iloop
	bbav=amax1(eqa,bbav/iloop)
	ebav=sqrt(ebav)/iloop
	if(.not.doplot) goto 81
	eav=eav/ego
	range=range/ego
	write (3,309) eav,range,bbav,ebav
81	if(.not. dolist) goto 12
	ln=ln+1
	bl(ln,iq)=bbav
	be(ln,iq)=ebav
	blen(ln,1)=eav-range
	blen(ln,2)=eav+range
12	e=e+del
	if(.not. doplot) goto 26
c
c	write out theory
1	if(.not. doplot) goto 26
	if (noth.eq.0) goto 30
	igr=jpart(iq)
c	temporary edit til vm1 and vm2 defined for kaon and proton beams
	if (igr.le.4) pll=plim2(1)
	if (igr.ge.5) pll=plim2(igr-4)
	ello=eplx-1.0
	note=50
	if (fixed(iq)) note=100
	if (nfs) note=200
	ddt=(ephx-eplx+2.)/(float(note)-1.)
	call wgtint(0.,ddt,moint,val1,wgt1)
	it=1
	ibm=(nobeam+1)/2
	call prepth(note,ddt,ello,pll,fixed(iq),igr)
c	if (igr.eq.5) call pioncs(note,ddt,ello,acpi)
	iline=0
	write(3,311) iline
	e=ello+0.5*ddt
	fudmt=fudm/ddt
	do 4 i=1,note
	thiq=th(i)*fudmt
	if (igr.eq.5) thiq=thiq+acpi(i)
	emma=e/ego
	if(thiq.gt.1.e-12) write(3,310) emma,thiq
4	e=e+ddt
c
	if (noth.le.1) goto 28
	do 5 j=1,noth
	if (j.eq.1.and.fixed(1).and.(.not.fixed(iq))) goto 5
	write(3,311) j
	e=ello+0.5*ddt
	thmax(j)=0.
	acth(j)=0.
	do 6 i=1,note
	thip=thy(i,j)*fudmt
	acth(j)=acth(j)+thip
	thmax(j)=amax1(thmax(j),thip)
	emma=e/ego
	if (th(i).lt.1.e-12.or.thip.lt.1.e-10) goto 95
	if (thy(i,j)/th(i).lt.0.999) write(3,310) emma,thip
	goto 6
95	if (i.eq.note) goto 6
	if (thy(i+1,j)*fudmt.lt.1.e-10.and.thy(i-1,j)*fudmt.lt.1.e-10)
&	goto 6
	if (th(i).lt.1.e-6) goto 6
	if (thip.gt.1.e-10) write(3,310) emma,thip
6	e=e+ddt
5	continue
	print 404,(acth(j),j=1,noth)
c
c	index to theory
28	if (iq.eq.1.and.iini2.gt.0) goto 34
	if (.not.fitalp(1)) goto 34
	alp=param(ialpha(1))
	alperr=perr(ialpha(1))
	salp=-alp
	if (iini2.eq.0) ccyy=ccyy-50.0
	if (alperr.gt.0.10) goto 77
c	if (alperr.gt.1.0.or.alp.ge.1.0.or.alp.le.-1.0) write(3,324) nsiz,
c&	ccl,ccyy,alp,alperr
c	if ((alp.lt.1.0.and.alp.gt.0.0).and.alperr.lt.1.0) write(3,327)
c&	nsiz,ccl,ccyy,alp,alperr
c	if ((alp.gt.-1.0.and.alp.lt.0.0).and.alperr.lt.1.0) write(3,328)
c&	nsiz,ccl,ccyy,salp,alperr
	goto 34
77	continue
c	if (alperr.gt.1.0.or.alp.ge.1.0.or.alp.le.-1.0) write(3,347) nsiz,
c&	ccl,ccyy,alp,alperr
c	if ((alp.lt.1.0.and.alp.gt.0.0).and.alperr.lt.1.0) write(3,348)
c&	nsiz,ccl,ccyy,alp,alperr
c	if ((alp.gt.-1.0.and.alp.lt.0.).and.alperr.lt.1.0) write(3,349)
c&	nsiz,ccl,ccyy,salp,alperr
34	continue
c	plot arrows
	xbeg=e1(iini)/ego
	xend=e2(iini)/ego
	if (iini1.gt.0.and.iq.eq.2) goto 36
	ipz=max0(ifix((e1(iini)-ello)/ddt)+1,1)
	ipy=max0(ifix((e2(iini)-ello)/ddt)+1,1)
	iea=1+ifix((e2(iini)-elow(iini))/ddi)
	iec=1+ifix((e1(iini)-elow(iini))/ddi)
	ybeg=0.
	yend=0.
	if (e1(iini).lt.eplx) goto 53
64	ybeg=th(ipz)*exp(1.0)*fudmt
	if (iec.ge.istrt) ybeg=amax1(ybeg,bnew(iec,iq)+enew(iec,iq))
	if (iec+1.ge.istrt) ybeg=amax1(ybeg,bnew(iec+1,iq)+enew(iec+1,iq))
	if (ybeg*expmul.gt.0.85*epa) ybeg=0.85*epa/expmul
	yend=ybeg*expmul
53	zbeg=th(ipy)*exp(1.0)*fudmt
	if (iea.le.noet) zbeg=amax1(zbeg,bnew(iea,iq)+enew(iea,iq))
	if (iea.lt.noet) zbeg=amax1(zbeg,bnew(iea+1,iq)+enew(iea+1,iq))
	if (iini1.eq.0.or.iini2.eq.0) goto 66
	ieb=1+ifix((e2(iini2)-elow(iini2))/de(iini2))
	if (ieb.le.noe(iini2)) zbeg=amax1(zbeg,bnew(ieb,2)+enew(ieb,2))
	if (ieb.lt.noe(iini2)) zbeg=amax1(zbeg,bnew(ieb+1,2)+
&	enew(ieb+1,2))
66	if (zbeg*expmul.gt.0.85*epa) zbeg=0.85*epa/expmul
	zend=zbeg*expmul
36	continue
c	write(3,331) xbeg,yend,xbeg,ybeg,xend,zend,xend,zbeg
30	continue
c	write(3,307)
26	continue
c
	if(.not. dolist) return
c
	j=0
	if ((lwbias(1).eq.0).and.(lwbias(2).eq.0)) goto 73
	if (ln.eq.0) goto 73
	if ((lwbias(2).eq.0).or.nfs) only1=.true.
cnono	call gutten(blen,bl,be,only1,ln,20,j)
	j=1
73	continue
cnono	call gutten(ben,b,ee,nfs,n,100,j)
c
c
	inpar=ipart
	nfsq=nfs
	thii=abs(thi)
	tloo=abs(tlo)
	avtt=abs(avt)
	if((lwbias(1).ne.0) .and. (lwbias(2).eq.0)) lwprob=.true.
cnono	call blycei(jpart)
	return
c
c	initialization
93	nplot=9
c	open files and rewind
	open(3,file='batchin',access='sequential',form='formatted')
	if (dolist) open(7,file='50OUT',access='sequential',form='formatted')
	if (dolist) open(9,file='52OUT',access='sequential',form='formatted')
	rewind(3)
	if (dolist) rewind(7)
	if (dolist) rewind(9)
c******write statements for first page commented out because a bug in
c******the plotting package does not allow plot pages without data
c	write(3,342)
	if (lwbias(1).eq.0.and.(.not.nfs)) ccxl=448.0
	labg=labely
	ccxp=ccxl+30.0
	if (noth.eq.0) goto 60
	do 7 j=1,noth
c	if (isymb(j).le.99) write(3,314) ccxp,labg,isymb(j),(nampl(i,j),
c&	i=1,3)
c	if (isymb(j).ge.100) write(3,340) ccxp,labg,isymb(j),(nampl(i,j),
c&	i=1,3)
7	labg=labg-40
c	write(3,315)
60	expmul=exp(50.0/183.3)
	write(3,316)
	nlab=11
	ntkx=5
	if (nfs) then
c	write(3,341)
	nono=7h$2x$u00
	else
	nono=3h$2x
	endif
	goto 94
c
	end
