* zcat ../output/u-000*.gz | awk 'BEGIN {e1=log(6.e2); e2=log(6.e4);} /^TR .* .* mu[+-] / {i=int(cos(3.141592654*$8/180)*100); j=int((log($(11))-e1)*100/(e2-e1)); a[i" "j]++} END {for(i=0;i<100;i++) for(j=0;j<100;j++) print a[i" "j]+0}' > ball.q2

* echo -e "\nexec muflux 0\nq" | paw | grep Res | awk '{print $4, $6, $7, $(14)}' > fit.f00

macro 1

do i=0.,0.99,0.01
exec a 1 2 std [1] 0 [i]
enddo

return


macro costh
exec init1
set XVAL 0.5

vec/cre lpar(1) R 1

vec/cre xhpar(5) R 0.00851164 0.124534 0.059761 2.32876 19.279
vec/cre xcpar(5) R 0.102573 -0.068287 0.958633 0.0407253 0.817285

appl comis quit
	real function f(x)
	vector lpar,xhpar,xcpar
	real aux, xs
	real h,c
	dimension par(5)
	i=lpar(1)

	if i.ge.3 then
		if x.lt.0.002 then
			par(1)=0.11137
			par(2)=0.
			par(3)=0.
			par(4)=0.
			par(5)=0.
		elseif x.lt.0.2 then
			par(1)=0.11148
			par(2)=-0.03427
			par(3)=5.2053
			par(4)=-14.197
			par(5)=16.138
		else
			par(1)=0.06714
			par(2)=0.71578
			par(3)=0.42377
			par(4)=-0.19634
			par(5)=-0.021145
		endif
	endif

	if i.eq.0 then
	xs=x
	elseif i.eq.1 then
	xs=sqrt((x*x+xcpar(1)*xcpar(1)+xcpar(2)*x**xcpar(3)+xcpar(4)
     *		*x**xcpar(5))/(1+xcpar(1)*xcpar(1)+xcpar(2)+xcpar(4)))
	elseif i.eq.2 then
		xs=sqrt(1-x*x)
		y=xs/(1+xhpar(1))
	aux=xhpar(5)*(1+xhpar(3)*xs**xhpar(4))/(1-(y*y))**xhpar(2)
		xs=xs/(1+aux/6.4E3)
		xs=sqrt(1-xs*xs)
	elseif i.ge.3 then
		sum=par(1)
		aux=1
		do n=2,5
		aux=aux*x
		sum=sum+par(n)*aux
		enddo
		xs=sum
	endif
	if i.eq.4 then
		aux=sqrt(x*x+1./400)
		xs=xs*(0.986+0.014/aux)
	endif
	f=xs
	end
quit

null 0 1 0 1

vec/inp lpar 0
set dmod 3
set hcol 7
fun/pl f 0 1 sl

vec/inp lpar 1
set dmod 1
set hcol 4
fun/pl f 0 1 sl

vec/inp lpar 2
set dmod 3
set hcol 3
fun/pl f 0 1 sl

vec/inp lpar 3
set dmod 2
set hcol 2
fun/pl f 0 1 sl

vec/inp lpar 4
set dmod 3
set hcol 2
fun/pl f 0 1 sl

null 0 1 0 1 sa

set TXCI 7
set PLCI 7
set dmod 3
dline 0.05 0.13 0.91 0.91
itx 0.15 0.90 no correction

set TXCI 4
set PLCI 4
set dmod 1
dline 0.05 0.13 0.86 0.86
itx 0.15 0.85 density ratio model

set TXCI 3
set PLCI 3
set dmod 3
dline 0.05 0.13 0.81 0.81
itx 0.15 0.80 local muon cos([q]) model

set TXCI 2
set PLCI 2
set dmod 2
dline 0.05 0.13 0.76 0.76
itx 0.15 0.75 Klimushin et al. result

set TXCI 2
set PLCI 2
set dmod 3
dline 0.05 0.13 0.71 0.71
itx 0.15 0.70 with sec[q] correction

set TXCI 1
atit 'cos([q])' 'cos([q^*!])'

exec plot1 costh
return


* ../ammc -p=gen/ExpCorr.java -c -r > out

macro mpa
exec init1

w=1
vec/cre lw(1) R [w]

appl comis quit
	function f(x)
	vector lw
	real x,y
	common /pawpar/ par(5)
	w=lw(1)
	if w.eq.1 then
		xs=sqrt(1-x*x)
		y=xs/(1+par(1))
		f=par(5)*(1+par(3)*xs**par(4))/(1-(y*y))**par(2)
	elseif w.eq.2 then
		f=sqrt((x*x+par(1)*par(1)+par(2)*x**par(3)
     *		+par(4)*x**par(5))/(1+par(1)*par(1)+par(2)+par(4)))
	elseif w.eq.3 then
		f=par(1)+par(2)*x**par(3)+par(4)*(1-x*x)**par(5)
	elseif w.eq.4 then
		f=par(1)+par(2)*x**par(3)+par(4)*(1-x*x)**par(5)
	endif
	end
quit

vec/cre chi2(1) R 0
appl comis quit
        subroutine chi2r
        COMMON/HCFITS/NCFITS,NPFITS,NFPAR,FITCHI,FITPAR(35),FITSIG(35),
     *          FITDER(35)
        vector chi2
        chi2(1)=sqrt(FITCHI)
        write(*,*) 'Chi2=', FITCHI
        end
quit

opt linx
opt liny
set XLAB 1.5
set XVAL 0


vec/rea x,h,c,o,l out4

sigma a=x
if [w].eq.1 then
sigma b=h
null 0 1 14 41
elseif [w].eq.2 then
sigma b=c
null 0 1 0 1
set dmod 3
set PLCI 7
dline 0 1 0 1
elseif [w].eq.3 then
sigma b=100/o
null 0 1 0 0.12
set XLAB 2.0
elseif [w].eq.4 then
sigma b=1000/l
null 0 1 0 70
endif

set dmod 1
set PLCI 7
vec/pl b%a ! sl


vec/rea x,h,c,o,l out

sigma a=x
if [w].eq.1 then
sigma b=h
elseif [w].eq.2 then
sigma b=c
set dmod 3
set PLCI 7
dline 0 1 0 1
elseif [w].eq.3 then
sigma b=100/o
set XLAB 2.0
elseif [w].eq.4 then
sigma b=1000/l
endif

set dmod 1
set PLCI 3
vec/pl b%a ! sl

vec/cre e(101) R 101*0.01

if [w].eq.1 then
vec/cre par(5) R 0.01 0.2 0 1 20
f1=2
f2=4
elseif [w].eq.2 then
vec/cre par(5) R 0.1 0 1 0 1
f1=3
f2=5
elseif [w].eq.3 then
vec/cre par(5) R 0.003 0.1 1 0 1
f1=3
f2=5
elseif [w].eq.4 then
vec/cre par(5) R 1.5 50 1 0 1
f1=3
f2=5
endif

lout=b(101)
if [w].eq.1 then
vec/inp par(5) [lout]
endif

appl comis quit
	subroutine sx
	vector par
	common /pawpar/ xpar(5)

	do i=1,5
		xpar(i)=par(i)
	enddo

	end
quit
appl comis quit
	call sx
	end
quit

sigma e=b/10000
set dmod 2
set PLCI 4
vec/fit a b e f se [f1] par
sigma par3=par

appl comis quit
        call chi2r
        end
quit
chi23=chi2(1)

if [w].eq.1 then
vec/inp par(3) 0.1
elseif [w].eq.2 then
vec/inp par(4) 0.01
elseif [w].eq.3 then
vec/inp par(4) 0.1
elseif [w].eq.4 then
vec/inp par(4) 0.1
endif

set dmod 3
set PLCI 2
vec/fit a b e f se [f2] par
sigma par5=par

appl comis quit
        call chi2r
        end
quit
chi25=chi2(1)

do i=1,5
pa=p//[i]
[pa]=par([i])
sa=s//[i]
qa=q//[i]

if [[pa]]<0 then
[sa]='-'
[qa]=-[[pa]]
else
[sa]='+'
[qa]=[[pa]]
endif

enddo

MESS [q1]*(1[s4][q4]*x**[p5])/(1-(x/(1[s2][q2]))**2)**[p3]

set hwid 5

if [w].eq.1 then
atit 'cos([q])' 'muon production altitude "M# km "N#'
elseif [w].eq.2 then
atit 'cos([q])' 'cos([q])^*!'
elseif [w].eq.3 then
atit 'cos([q])' '1/mass overburden "M# 1/mwe "N#'
elseif [w].eq.4 then
atit 'cos([q])' '1/muon track length "M# 10^-3! 1/km "N#'
endif

null 0 1 0 1 sa

if [w].eq.1 then
set chhe 0.4
elseif [w].eq.2 then
set chhe 0.3
elseif [w].eq.3 then
set chhe 0.4
elseif [w].eq.4 then
set chhe 0.4
endif

set TXCI 3
set PLCI 3
set dmod 1
dline 0.08 0.16 0.91 0.91
itx 0.18 0.90 model with X?0!=114.8 g/cm^2!, h?0!=19.28 m

set TXCI 4
set PLCI 4
set dmod 2
dline 0.08 0.16 0.86 0.86
chi23=[chi23]/100
chi2a=$FORMAT([chi23] ,f5.2)
itx 0.18 0.85 '3-parameter fit to the model (av. dev. '//[chi2a]//'%)'

set TXCI 2
set PLCI 2
set dmod 3
dline 0.08 0.16 0.81 0.81
chi25=[chi25]/100
chi2a=$FORMAT([chi25] ,f5.2)
itx 0.18 0.80 '5-parameter fit to the model (av. dev. '//[chi2a]//'%)'

set TXCI 7
set PLCI 7
set dmod 1
dline 0.08 0.16 0.76 0.76
itx 0.18 0.75 model for the South Pole July atmosphere

MESS

set TXCI 1
set PLCI 1
set chhe 0.4
MESS

MESS 3-parameter fit:
vec/pri par3
MESS

MESS 5-parameter fit:
vec/pri par5
MESS

MESS use for CORSIKA fit:
p1=par5(1)
p2=par5(2)
p3=par5(3)
p4=par5(4)
p5=par5(5)
MESS par[w] [p1] [p2] [p3] [p4] [p5]

exec plot1 'mpa'//[w]
return


* zcat ../output/u-000* | awk '/^TR .* .* mu[+-] / {i=int(cos(3.141592654*$8/180)*100); a[i]+=$(10); n[i]++} END {for(i=0;i<100;i++) print (i+0.5)/100, a[i]/n[i]}' > zl

macro zl
exec init1

set MTYP 7
SET MSCF 0

vec/rea x,h,c,o,l out
sigma l=1000/l
vec/rea q,p zl
sigma p=1000000/p
set XLAB 1.5
set XVAL 0

null 0 1 0 70
set PLCI 3
set dmod 1
vec/pl l%x ! sl
set PLCI 4
set dmod 2
vec/pl p%q ! sl

atit 'cos([q])' '1/muon track length "M# 10^-3! 1/km "N#'

null 0 1 0 1 sa
set chhe 0.4

set TXCI 3
set PLCI 3
set dmod 1
dline 0.08 0.16 0.91 0.91
itx 0.18 0.90 model with X?0!=114.8 g/cm^2!, h?0!=19.28 m

set TXCI 4
set PLCI 4
set dmod 2
dline 0.08 0.16 0.86 0.86
itx 0.18 0.85 average muon track length from CORSIKA

exec plot1 zl
return


macro atm
exec init1

set MTYP 7
SET MSCF 0

vec/rea gh,go,gh2,gr,gr2,gh0 atm1
vec/rea ch,co,ch2,cr,cr2,ch0 atm2
vec/rea sh,so,sh2,sr,sr2,sh0 atm3

zone 1 2

opt logy
null 0 100 1.e-4 1.e4

set dmod 1
set PLCI 2
vec/pl go%gh ! sl

set dmod 2
set PLCI 4
vec/pl co%ch ! sl

set dmod 3
set PLCI 3
vec/pl so%sh ! sl

set dmod 1
set PLCI 2
vec/pl gr%gh ! sl

set dmod 2
set PLCI 4
vec/pl cr%ch ! sl

set dmod 3
set PLCI 3
vec/pl sr%sh ! sl

opt liny
null 0 1 0 1 sa
set dmod 1
set PLCI 2
set TXCI 2
dline 0.2 0.32 0.93 0.93
itx 0.35 0.9 Shibata (from Gaisser, 1990)
set dmod 2
set PLCI 4
set TXCI 4
dline 0.2 0.32 0.83 0.83
itx 0.35 0.8 Standard US (as in CORSIKA)
set dmod 3
set PLCI 3
set TXCI 3
dline 0.2 0.32 0.73 0.73
itx 0.35 0.7 South Pole July 1 (CORSIKA)

set TXCI 1
atit '' 'X "M# g/cm^2! "N# and dX/dh!'

opt liny
null 0 100 0. 20

set dmod 1
set PLCI 2
vec/pl gh0%gh ! sl

set dmod 2
set PLCI 4
vec/pl ch0%ch ! sl

set dmod 3
set PLCI 3
vec/pl sh0%sh ! sl

atit 'height "M# km "N#' 'h?0! "M# km "N#'

exec plot1 atm
return


macro a
exec init1

presm=1
presf=3

* partl: mu=1 nu_mu=2 nu_e=3
partl=1
* atmsf: std_us_0=0 spl_oct=1 std_us_2834=2 spl_jul_2834=3 spl_jul_-0.115548=4
atmsf=0

thx=0.1
tjx=0.

presm=[4]
presf=[5]
thx=[6]
*tjx=[6]
MESS [presm] [presf] [tjx] [thx] [partl] [atmsf]

th1=1-[thx]
th2=[thx]*100
th3=[th1]*100

tj1=1-[tjx]
tj2=[tjx]*100
tj3=[tj1]*100

* fluxsum=0.361788 * ws
fluxsum=0.323626

if [1].eq.' ' then
	mflag=1
else
	mflag=[1]
endif

if [2].eq.' ' then
	lflag=1
else
	lflag=[2]
endif

if [3].eq.' ' then
	mxll='std'
else
	mxll=[3]
endif

if [mflag].eq.1 then
	model=qgsjet
	corfact=1
	showers=100
	end='q'
endif

showers=[showers]*100*1.e6

if [lflag].eq.1 then
	emin=1.e3
	emax=1.e5
elseif [lflag].eq.2 then
	emin=6.e2
	emax=6.e4
	end=[end]//'2'
endif

MESS model=[model] (mflag=[mflag]) lflag=[lflag] mxll=[mxll] file='b'//[mxll]//'.'//[end]

lmin=$SIGMA(log([emin])/log(10))
lmax=$SIGMA(log([emax])/log(10))
ldif=[lmax]-[lmin]
lsum=[lmax]+[lmin]
vec/cre lpar(2) R [lsum] [presm]

fmin=1
fmax=1.e6

scale=$SIGMA(10000/log(10)/[ldif])
scale=[scale]*([fluxsum]/[showers])*[corfact]*1.e-4

vec/cre a(100,100) R
vec/cre b(100,100) R
vec/cre c([th3],100) R
vec/cre cj([tj3],100) R
vec/rea a 'b'//[mxll]//'.'//[end]

do i=1,100
do j=1,100
aux=a([i],[j])
ii=101-[i]
vec/inp b([j],[ii]) [aux]
if [j].gt.[th2] then
jj=[j]-[th2]
vec/inp c([jj],[ii]) [aux]
endif
if [j].gt.[tj2] then
jj=[j]-[tj2]
vec/inp cj([jj],[ii]) [aux]
endif
enddo
enddo
sigma a=b

2d 100 ' ' 100 0 1 100 [lmin] [lmax]
hi/put/cont 100 a
sigma e=sqrt(a)
hi/put/err 100 e

2d 300 ' ' [th3] [thx] 1 100 [lmin] [lmax]
hi/put/cont 300 c
sigma e=sqrt(c)
*sigma e=a/100
hi/put/err 300 e

2d 400 ' ' [tj3] [tjx] 1 100 [lmin] [lmax]
hi/put/cont 400 cj
sigma e=sqrt(cj)
*sigma e=a/100
hi/put/err 400 e

if [atmsf].eq.0 then
* CORSIKA standard US atmosphere 0 m elevation
vec/cre xhpar(5) R 0.00851164 0.124534 0.059761 2.32876 19.279
vec/cre xcpar(5) R 0.102573 -0.068287 0.958633 0.0407253 0.817285
vec/cre xopar(5) R -0.017326 0.114236 1.15043 0.0200854 1.16714
vec/cre xlpar(5) R 1.3144 50.2813 1.33545 0.252313 41.0344
elseif [atmsf].eq.1 then
* CORSIKA SPl OCT01 atmosphere 2.834 m elevation
vec/cre xhpar(5) R 0.00762702 0.143502 0.0806641 2.32839 14.9332
vec/cre xcpar(5) R 0.0869606 -0.162557 1.17013 0.0852275 1.00028
vec/cre xopar(5) R -0.0255075 0.168142 1.1455 0.0292488 1.18082
vec/cre xlpar(5) R 0.957459 65.6647 1.37996 0.743004 8.70671
elseif [atmsf].eq.2 then
* CORSIKA standard US atmosphere 2.834 m elevation
vec/cre xhpar(5) R 0.00824021 0.138291 0.0798717 2.3316 16.4457
vec/cre xcpar(5) R 0.0987167 -0.0667872 0.960375 0.0387649 0.811397
vec/cre xopar(5) R -0.024379 0.161448 1.148 0.0281572 1.17219
vec/cre xlpar(5) R 1.0135 59.4791 1.37163 0.629313 9.50148
elseif [atmsf].eq.3 then
* CORSIKA SPl JUL01 atmosphere 2.834 m elevation
vec/cre xhpar(5) R 0.00722581 0.136347 0.083569 2.3381 14.1322
vec/cre xcpar(5) R 0.0937706 -0.0222859 0.577757 0.0140503 0.458109
vec/cre xopar(5) R -0.0263716 0.171436 1.14647 0.0301484 1.18237
vec/cre xlpar(5) R 0.884797 69.5992 1.37278 0.878282 6.0362
elseif [atmsf].eq.4 then
* CORSIKA SPl JUL01 atmosphere 0 m elevation
vec/cre xhpar(5) R 0.00750474 0.120382 0.0593859 2.32484 17.0799
vec/cre xcpar(5) R 0.098064 -0.0191912 0.570608 0.0114115 0.435399
vec/cre xopar(5) R -0.0175821 0.114603 1.14833 0.0201985 1.17589
vec/cre xlpar(5) R 1.45034 56.799 1.32806 0.214639 57.6376
endif

if [partl].eq.1 then
	vec/cre xvpar(3) R 115. 850. 0.054
	a014=0.14
elseif [partl].eq.2 then
	a1=121*1.1/6
	a2=897*1.1/1.44
	a3=0.213
	a014=2.85e-2
	vec/cre xvpar(3) R [a1] [a2] [a3]
elseif [partl].eq.3 then
	vec/cre xvpar(3) R 0 0 0
	a014=2.4e-3
endif


vec/cre vc(8) R 1 0.6 0.4 0.3 0.2 0.1 0.05 0
vec/cre va(8) R -1 -0.355 -0.687 -0.619 -0.384 -0.095 0 0.083
vec/cre vb(8) R 0 -0.23 -0.01 -0.007 -0.09 -0.165 -0.186 -0.215

appl comis quit
	real function gaisser(x,y)
	implicit none
	vector lpar,xcpar,xopar,xlpar,xvpar,va,vb,vc
	real x,y,aux, ei, ef, aux1, aux2, auf1, auf2
	real c,o,l,a,b,as,bs,p,de,zfi,zff,f1,f2,f3,f4,vf
	integer flag, i
	common /pawpar/ par(5)
	real gpar
	dimension gpar(3)
	flag=lpar(2)
	vf=0

	ef=10**(lpar(1)-y)
	do i=1,3
		gpar(i)=(exp(par(i+2))+10)/(10*exp(par(i+2))+1)
	enddo

	if flag.ge.0 then
	c=x
	o=0
	l=0
	ei=ef
	de=1
	p=1
	endif

	if flag.ge.1 then
	aux=xcpar(1)*xcpar(1)
	c=sqrt((x*x+(aux+xcpar(2)*x**xcpar(3)+xcpar(4)*x**xcpar(5))*
	gpar(1))/(1+(aux+xcpar(2)+xcpar(4))*gpar(1)))
	endif

	if flag.ge.2 then
	o=1/(xopar(1)+xopar(2)*x**xopar(3)+xopar(4)*(1-x*x)**xopar(5))

	o=(o*gpar(2)-1.148)
	if o.lt.0 then
		o=0
	endif
	a=0.262
	b=0.350e-3
	ei=((a+b*ef)*exp(b*o)-a)/b

	as=-0.487
	bs=8.766e-3

	aux1=(bs*a-as*b)
	aux1=2*aux1*aux1
	aux2=2*b*b*b
	zfi=(bs*ei*b*(-2*bs*a+4*as*b+bs*ei*b)+aux1*log(a+ei*b))/aux2
	zff=(bs*ef*b*(-2*bs*a+4*as*b+bs*ef*b)+aux1*log(a+ef*b))/aux2

	aux=1.1*c
	auf1=aux/xvpar(1)
	auf2=aux/xvpar(2)
	aux1=1/(1+auf1*ei)
	aux2=1/(1+auf2*ei)
	aux=aux1+xvpar(3)*aux2

	f1=ei**(-par(2)-1)*aux
	f2=ei**(-par(2)-1)*(-(par(2)+1)*aux/ei
		-auf1*aux1*aux1-xvpar(3)*auf2*aux2*aux2)
	f3=(-par(2)-1)*f2/ei+ei**(-par(2)-1)*
		((par(2)+1)*aux/(ei*ei)+((par(2)+1)/ei)*
		(auf1*aux1*aux1+xvpar(3)*auf2*aux2*aux2)
		+2*auf1*auf1*aux1*aux1*aux1
		+2*xvpar(3)*auf2*auf2*aux2*aux2*aux2)
	aux=f2/f1
	f4=f3/f1-aux*aux

	aux1=as+bs*ef
	aux2=as+bs*ei
	de=exp(b*o)-(f2/(2*f1))*(aux1*aux1-aux2*aux2)/(a+b*ef)-
		(f4/2)*(zff-zfi)

	ei=ei-(f2/(2*f1))*(zff-zfi)
	endif

	if flag.ge.3 then
	l=1E3/(xlpar(1)+xlpar(2)*x**xlpar(3)+xlpar(4)*(1-x*x)**xlpar(5))
	p=l*0.105658357/(2.19703e-6*299792.458*ei)
	p=exp(-p*gpar(3))
	endif

	if xvpar(1).gt.0 then
		aux=xvpar(3)/(1+1.1*ei*c/xvpar(2))
		aux=aux+1/(1+1.1*ei*c/xvpar(1))
		gaisser=p*de*par(1)*ei**(-par(2))*aux*ef/ei
	else
		if vf.eq.0 then
		if x.lt.0.3 then
			f1=0.11-2.4*x
			f2=-0.22+0.69*x
		else
			f1=-0.46-0.54*x
			f2=-0.01+0.01*x
		endif
		else
		if x.gt.vc(2) then
			i=2
		elseif x.gt.vc(3) then
			i=3
		elseif x.gt.vc(4) then
			i=4
		elseif x.gt.vc(5) then
			i=5
		elseif x.gt.vc(6) then
			i=6
		elseif x.gt.vc(7) then
			i=7
		else
			i=8
		endif
		f1=va(i)+(va(i-1)-va(i))*(x-vc(i))/(vc(i-1)-vc(i))
		f2=vb(i)+(vb(i-1)-vb(i))*(x-vc(i))/(vc(i-1)-vc(i))
		endif

		aux=11.4*ei**(f1+f2*log(ei)/log(10.))
		aux=aux/(1+1.21*ei*c/121)
		aux=aux+0.05/(1+1.5*ei*c/897)
		aux=aux+0.185/(1+1.5*ei*c/194)
		gaisser=p*de*par(1)*ei**(-par(2))*aux*ef/ei

	endif
	end
quit

vec/cre xpar(5) R 1.e10 1.72 1 1 1

sigma xap=log((10-xpar)/(10*xpar-1))
do i=3,5
ap=xap([i])
vec/inp xpar([i]) [ap]
enddo

appl comis quit
	subroutine parout
	common /pawpar/ par(5)
	vector xpar
	do i=1,5
		par(i)=xpar(i)
	enddo
	end
	call parout
	end
quit

vec/cre chi2(1) R 0
appl comis quit
        subroutine chi2r
        COMMON/HCFITS/NCFITS,NPFITS,NFPAR,FITCHI,FITPAR(35),FITSIG(35),
     *          FITDER(35)
        vector chi2
        chi2(1)=FITCHI
        write(*,*) 'Chi2=', FITCHI
        end
quit

hi/fit 300 gaisser s 1 xpar
hi/fit 300 gaisser se 2 xpar

if [presf].ge.1 then
presa=[presf]+2
hi/fit 300 gaisser se [presa] xpar
endif

hi/fit 400 gaisser s 1 xpar

appl comis quit
        call chi2r
        end
quit
chi2a=chi2(1)

hi/fit 300 gaisser s 1 xpar

appl comis quit
        call chi2r
        end
quit
chi2=chi2(1)

opt logz
fun2 200 gaisser 20 0 1 20 [lmin] [lmax] ' '

call hminim(100,[fmin])
call hmaxim(100,[fmax])
call hminim(200,[fmin])
call hmaxim(200,[fmax])

set ndvx 10505
set ndvy -1
set XVAL 999

2d 900 ' ' 1 0.001 1 1 [lmin] [lmax]
gmin=[fmin]*1.001
call hminim(900,[gmin])
call hmaxim(900,[fmax])

set hcol 8
hi/pl 900 'surf fb'

set hcol 4
hi/pl 100 'surf4 fb bb as'

set hcol 1
hi/pl 200 'surf fb bb as'

appl comis quit
	subroutine parout
	common /pawpar/ par(6)
	write(*,*) par
	end
	call parout
	end
quit

set TXCI 1
set XLAB 1.0
set YLAB 1.0

atit 'cos([q])' 'E "M#GeV"N#' '' 312

null 0 1 0 1 sa
set lwid 2

jxi=0
jxf=0.365
jyi=0.279
jyf=0

do jf=2,5
do jg=0,9

je=$SIGMA([jf]+log([jg])/log(10))
if ([je].ge.[lmin]).and.([je].le.[lmax]) then

jxx=[jxi]+([jxf]-[jxi])*([je]-[lmin])/([lmax]-[lmin])
jyy=[jyi]+([jyf]-[jyi])*([je]-[lmin])/([lmax]-[lmin])

ixx=[jxx]
if [jg].eq.0 then
iyy=[jyy]-0.02
else
iyy=[jyy]-0.01
endif

dline [ixx] [jxx] [iyy] [jyy]
endif

enddo
enddo

if [lflag].eq.1 then
	itx -0.04  0.18 10^3!
	itx  0.13  0.05 10^4!
	itx  0.30 -0.08 10^5!
elseif [lflag].eq.2 then
	itx 0.00 0.15 10^3!
	itx 0.17 0.02 10^4!
endif

xi=0.23
xd=0.127
yi=-0.10
yd=0.033

if [lflag].eq.1 then
	gg=0.2
	xi=[xi]+[xd]
	yi=[yi]+[yd]
elseif [lflag].eq.2 then
	gg=0.0
endif

do g=[gg],1.0,0.2
xi=[xi]+[xd]
yi=[yi]+[yd]
itx [xi] [yi] [g]
enddo

A=xpar(1)
A=[A]*[scale]/[a014]
g=xpar(2)
g=[g]+1

vec/pri xpar
MESS [A] [g]


sigma oupar=(exp(xpar)+10)/(10*exp(xpar)+1)
p1=oupar(3)
p2=oupar(4)
p3=oupar(5)
MESS
MESS Result: [presm] [presf] [thx] [tjx]  [A] [g] [p1] [p2] [p3]  [chi2] [chi2a]

itx 0.80 1.00 A
itx 0.81 0.94 '[g]'
itx 0.80 0.88 [h]^2!
itx 0.85 1.00 = $FORMAT([A],f5.3)
itx 0.85 0.94 = $FORMAT([g],f5.3)
itx 0.85 0.88 = $FORMAT([chi2] ,f5.2)

itx 0.00 1.00 [F]?0! = [a014] cm^-2!s^-1!sr^-1! [\327] A
set chhe 0.4
itx 0.00 0.94 '[h]^2!='//$FORMAT([chi2a] ,f5.2)//' (cos[q]=0-1)'
itx 0.00 0.88 [model]

appl comis quit
	real function gs(x)
	vector lpar
	real eng,aux
	eng=lpar(1)-4
	aux=gaisser(1.0, eng)
	gs=gaisser(x, eng)/aux
	end
quit

set ndvx
set ndvy
set XVAL
set XLAB 0
opt linx
opt liny
* null 0 1 0 13
set dmod 1
* fun/pl gs 0 1 s
* atit 'cos([q])' 'I(cos([q])) / I(cos([q])=1)'

vec/cre x(101) R
vec/cre y(101) R
do i=1,101
k=([i]-1)/100
vec/inp x([i]) [k]
ff=gs([k])
vec/inp y([i]) $CALL([ff])
enddo
vec/wri x,y gsx

exec plot1 'a'//[mflag]//[lflag]//[mxll]//[presm]//[presf]//[thx]
return


macro ab
exec init1
set XVAL 0

vec/cre c(7) R 1 0.4 0.3 0.2 0.1 0.05 0
vec/cre a(7) R -1 -0.687 -0.619 -0.384 -0.095 0 0.083
vec/cre b(7) R 0 -0.01 -0.007 -0.09 -0.165 -0.186 -0.215
vec/cre e(7) R 1 1 1 1 1 1 1

vec/cre c(8) R 1 0.6 0.4 0.3 0.2 0.1 0.05 0
vec/cre a(8) R -1 -0.355 -0.687 -0.619 -0.384 -0.095 0 0.083
vec/cre b(8) R 0 -0.23 -0.01 -0.007 -0.09 -0.165 -0.186 -0.215
vec/cre e(8) R 1 1 1 1 1 1 1 1

vec/cre c1(5) R 0.3 0.2 0.1 0.05 0
vec/cre a1(5) R -0.619 -0.384 -0.095 0 0.083
vec/cre b1(5) R -0.007 -0.09 -0.165 -0.186 -0.215
vec/cre e1(5) R 1 1 1 1 1

vec/cre c2(2) R 1 0.3
vec/cre a2(2) R -1 -0.619
vec/cre b2(2) R 0 -0.007
vec/cre e2(2) R 1 1

sigma e=e/100

zone 1 2
null 0 1 -1.1 0.1
set dmod 1
set PLCI 4
vec/pl a%c ! sl
set dmod 2
set PLCI 2
vec/fit c1 a1 e P1 s
vec/fit c2 a2 e P1 s
atit ! 'a'

null 0 1 -0.25 0.05
set dmod 1
set PLCI 4
vec/pl b%c ! sl
set dmod 2
set PLCI 2
vec/fit c1 b1 e P1 s
vec/fit c2 b2 e P1 s
atit 'cos[q]' 'b'

sigma e=a+b*log(600)/log(10)

exec plot1 ab
return


macro b
exec init1

set XLAB 1.5
set XVAL 0

null 0 1 0 13
set dmod 1
set MTYP 7
SET MSCF 0

vec/rea a,b gsnc
vec/inp b(101) 1
set dmod 3
set PLCI 3
set PMCI 3
vec/pl b(8:101)%a(8:101) ! sl

vec/rea a,b gsc1
vec/inp b(101) 1
set dmod 2
set PLCI 7
set PMCI 7
vec/pl b%a ! sl

vec/rea a,b gsc3
vec/inp b(101) 1
set dmod 1
set PLCI 4
set PMCI 4
vec/pl b%a ! sl

vec/rea a,b gsc4
vec/inp b(101) 1
set dmod 2
set PLCI 4
set PMCI 4
vec/pl b%a ! sl

set lwid 5
atit 'cos([q])' 'I(cos([q])) / I(cos([q])=1) at 10 TeV'


null 0 1 0 1 sa
set chhe 0.4

set TXCI 3
set PLCI 3
set dmod 3
dline 0.20 0.28 0.91 0.91
itx 0.30 0.90 no correction

set TXCI 7
set PLCI 7
set dmod 2
dline 0.20 0.28 0.86 0.86
itx 0.30 0.85 cos([q])-correction only

set TXCI 4
set PLCI 4
set dmod 1
dline 0.20 0.28 0.81 0.81
itx 0.30 0.80 fit with the corrected formula

set TXCI 4
set PLCI 4
set dmod 2
dline 0.20 0.28 0.76 0.76
itx 0.30 0.75 for the South Pole July atmosphere

exec plot1 b
return

* zcat ../output/u-000* | awk '/^TR .* .* mu[+-] / {i=int(cos(3.141592654*$8/180)*1000); a[i]++} END {for(i=0;i<1000;i++) print (i+0.5)/1000, a[i]+0}' > zenith.1
* zcat ../output/u-000* | awk '/^TR .* .* mu[+-] / {if($(11)>10000){ i=int(cos(3.141592654*$8/180)*1000); a[i]++}} END {for(i=0;i<1000;i++) print (i+0.5)/1000, a[i]+0}' > zenith.02
* cat zenith.02 | awk '{a+=$1; b+=$2} NR%5==0 { print a/5, b; a=0; b=0 }' > zenith.2

macro z
exec init1

flag=1

fluxsum=0.323626
corfact=1
showers=100
showers=[showers]*100*1.e6

scale=([fluxsum]/[showers])*[corfact]
MESS [scale]

if [flag].eq.1 then
	max=0.45e-5*0.14
	file='.1'
elseif [flag].eq.2 then
	file='.2'
	max=0.50e-8*0.14
endif
mxx=10*[max]/[scale]

if [flag].eq.1 then
	mnx=[max]*1.e6
elseif [flag].eq.2 then
	mnx=[max]*1.e9
endif

if [flag].eq.2 then
	mxx=[mxx]*5
endif

set XVAL 0
null 0 1 0 [mnx]

null 0 1 0 [mxx] sa
vec/rea a,b zenith//[file]
vec/pl b%a ! sl

set dmod 1
set PLCI 3
set lwid 10
set PMCI 0
null 0 1 0 [max] sa
vec/rea x,z zfit//[file]
vec/pl z%x ! sl
set lwid 5

set dmod 2
set PLCI 4
set PMCI 0
null 0 1 0 [max] sa
vec/rea x,z zf17//[file]
vec/pl z%x ! sl

set dmod 3
set PLCI 2
set PMCI 0
null 0 1 0 [max] sa
vec/rea x,z zf1p//[file]
vec/pl z%x ! sl

set PLCI 5
set dmod 2
a1=[max]*0.7
dline 0.3 0.3 0 [a1]

null 0 1 0 1 sa
set chhe 0.4

set TXCI 1
set PLCI 1
set dmod 1
dline 0.20 0.28 0.91 0.91
itx 0.30 0.90 CORSIKA muon integrated flux

set TXCI 3
set PLCI 3
set dmod 1
dline 0.20 0.28 0.86 0.86
itx 0.30 0.85 fit with the corrected formula

set TXCI 4
set PLCI 4
set dmod 2
dline 0.20 0.28 0.81 0.81
itx 0.30 0.80 fit with the original formula

set TXCI 2
set PLCI 2
set dmod 3
dline 0.20 0.28 0.76 0.76
itx 0.30 0.75 fit with cos([q])-correction only

set TXCI 1
if [flag].eq.1 then
	atit 'cos([q])' 'Muon flux above 600 GeV "M# 10^-6! cm^-1! sr^-1! s^-1! "N#'
elseif [flag].eq.2 then
	atit 'cos([q])' 'Muon flux above 10 TeV "M# 10^-9! cm^-1! sr^-1! s^-1! "N#'
endif

exec plot1 'z'//[flag]
return


macro fit
exec init1
set XVAL 0

flag=1

if [flag].eq.1 then
gmin=2.65
gmax=2.75
amin=0.50
amax=0.74
cmin=1
cmax=200
hmin=1
hmax=2
vec/rea x,a,g,c fit.f30
vec/rea x2,a2,g2,c2 fit.f00
vec/rea x3,a3,g3,c3 fit.f10
vec/rea x4,a4,g4,c4 fit.f20
elseif [flag].eq.2 then
gmin=2.70
gmax=2.73
amin=0.67
amax=0.74
cmin=1
cmax=20
hmin=1
hmax=1.2
vec/rea x,a,g,c fit.f30
vec/rea x2,a2,g2,c2 fit.h11
vec/rea x3,a3,g3,c3 fit.ha4
vec/rea x4,a4,g4,c4 fit.hmf
endif

set MTYP 7
SET MSCF 0

zone 1 3
null 0 1 [gmin] [gmax]
 set PLCI 4
set dmod 1
vec/pl g%x ! sl
g5=g(51)
set lwid 1
dline 0 1 [g5] [g5]
set lwid 5
 set PLCI 3
set dmod 4
vec/pl g2%x2 ! sl
g5=g2(51)
set lwid 1
dline 0 1 [g5] [g5]
set lwid 5
 set PLCI 7
set dmod 3
vec/pl g3%x3 ! sl
g5=g3(51)
set lwid 1
dline 0 1 [g5] [g5]
set lwid 5
 set PLCI 6
set dmod 2
vec/pl g4%x4 ! sl
g5=g4(51)
set lwid 1
dline 0 1 [g5] [g5]
set lwid 5
atit ! '[g]'

null 0 1 0 1 sa

if [flag].eq.1 then

set PLCI 3
set TXCI 3
set dmod 4
dline 0.02 0.13 1.23 1.23
itx 0.15 1.2 original formula

set PLCI 7
set TXCI 7
set dmod 3
dline 0.02 0.13 1.08 1.08
itx 0.15 1.05 +cos([q])-correction

set PLCI 6
set TXCI 6
set dmod 2
dline 0.57 0.68 1.23 1.23
itx 0.70 1.20 +dE/dx losses

set PLCI 4
set TXCI 4
set dmod 1
dline 0.57 0.68 1.08 1.08
itx 0.70 1.05 +muon decay

elseif [flag].eq.2 then

set PLCI 3
set TXCI 3
set dmod 4
dline 0.02 0.13 1.23 1.23
itx 0.15 1.2 fit in 1-100 TeV

set PLCI 7
set TXCI 7
set dmod 3
dline 0.02 0.13 1.08 1.08
itx 0.15 1.05 with SP Jul atm

set PLCI 6
set TXCI 6
set dmod 2
dline 0.52 0.63 1.23 1.23
itx 0.65 1.20 w/magnetic field

set PLCI 4
set TXCI 4
set dmod 1
dline 0.52 0.63 1.08 1.08
itx 0.65 1.05 corrected formula

endif

set TXCI 1

null 0 1 [amin] [amax]
 set PLCI 4
set dmod 1
vec/pl a%x ! sl
a5=a(51)
set lwid 1
dline 0 1 [a5] [a5]
set lwid 5
 set PLCI 3
set dmod 4
vec/pl a2%x2 ! sl
a5=a2(51)
set lwid 1
dline 0 1 [a5] [a5]
set lwid 5
 set PLCI 7
set dmod 3
vec/pl a3%x3 ! sl
a5=a3(51)
set lwid 1
dline 0 1 [a5] [a5]
set lwid 5
 set PLCI 6
set dmod 2
vec/pl a4%x4 ! sl
a5=a4(51)
set lwid 1
dline 0 1 [a5] [a5]
set lwid 5
atit ! 'A'

set XVAL 0.7
opt logy
null 0 1 [cmin] [cmax]
set PLCI 4
set dmod 1
vec/pl c%x ! sl
c5=c(51)
set lwid 1
dline 0.3 1 [c5] [c5]
set lwid 5
null 0 1 [hmin] [hmax] sa
set PLCI 4
vec/pl c(31:100)%x(31:100) ! sl
c5=c(51)
set lwid 1
dline 0.3 1 [c5] [c5]
set lwid 5

null 0 1 [cmin] [cmax] sa
set PLCI 3
set dmod 4
vec/pl c2%x2 ! sl
c5=c2(51)
set lwid 1
dline 0.3 1 [c5] [c5]
set lwid 5
null 0 1 [hmin] [hmax] sa
set PLCI 3
vec/pl c2(31:100)%x2(31:100) ! sl
c5=c2(51)
set lwid 1
dline 0.3 1 [c5] [c5]
set lwid 5

null 0 1 [cmin] [cmax] sa
set PLCI 7
set dmod 3
vec/pl c3%x3 ! sl
c5=c3(51)
set lwid 1
dline 0.3 1 [c5] [c5]
set lwid 5
null 0 1 [hmin] [hmax] sa
set PLCI 7
vec/pl c3(31:100)%x3(31:100) ! sl
c5=c3(51)
set lwid 1
dline 0.3 1 [c5] [c5]
set lwid 5

null 0 1 [cmin] [cmax] sa
set PLCI 6
set dmod 2
vec/pl c4%x4 ! sl
c5=c4(51)
set lwid 1
dline 0.3 1 [c5] [c5]
set lwid 5
null 0 1 [hmin] [hmax] sa
set PLCI 6
vec/pl c4(31:100)%x4(31:100) ! sl
c5=c4(51)
set lwid 1
dline 0.3 1 [c5] [c5]
set lwid 5
atit 'lower cos([q]) bound in fits' '[h]^2!'

opt liny
null 0 1 0 1 sa
itx 1.01 -0.05 [hmin]
itx 1.01 0.95 [hmax]

exec plot1 'fit'//[flag]
return


macro loss
exec init1


vec/cre par1(3) R

appl comis quit
        function f(x)
        common /pawpar/ par(3)
	f=(par(1)+par(2)*x)*(1-exp(-par(2)*par(3)))/(par(2)*par(3))
        end

	function zf(x)
	vector par1
        common /pawpar/ par(3)
	as=par(1)
	bs=par(2)
	a=par1(1)
	b=par1(2)
	zf=(bs*x*b*(-2*bs*a+4*as*b+bs*x*b)+
		2*(bs*a-as*b)**2*log(a+x*b))/(2*b**3)
	end

        function z(x)
	vector par1
        common /pawpar/ par(3)
	as=par(1)
	bs=par(2)
	a=par1(1)
	b=par1(2)
	o=par(3)
	ef=((a+b*x)*exp(-b*o)-a)/b
	z=sqrt((zf(x)-zf(ef))/o)
        end
quit

vec/cre xpar(3) R 1 1.e-3 10
appl comis quit
        subroutine parout
        common /pawpar/ par(3)
        vector xpar
        do i=1,3
                par(i)=xpar(i)
        enddo
        end
        call parout
        end
quit



opt grid
opt stat
opt logy
opt logx
null 100 1e7 1e-1 1e5

vec/rea a,b,c loss.1 ! 'OC'
sigma b=b/10
sigma c=c/sqrt(10)
num=$VDIM(a,1)

nnx=9
max=[num]-11
MESS BEGIN AT [nnx] END AT [max]

vec/cre er([num])
sigma er=b*0.001
vec/fit a([nnx]:[max]) b([nnx]:[max]) er([nnx]:[max]) f se0 2 xpar

appl comis quit
real function plot(x)
	common /pawpar/ par(2)
	plot=par(1)+par(2)*x
end
quit

set PLCI 6
set lwid 6
set dmod 1

annx=$SIGMA(a([nnx]))
amax=$SIGMA(a([max]))
b1nnx=$SIGMA(b([nnx])*3.)
b2nnx=$SIGMA(b([nnx])*.3)
b1max=$SIGMA(b([max])*3.)
b2max=$SIGMA(b([max])*.3)
dline [annx] [annx] [b1nnx] [b2nnx]
dline [amax] [amax] [b1max] [b2max]

set PLCI 7
set lwid 3
set dmod 1

f/pl f 1e-1 1e11 s
sigma par1=xpar

*fun/pl 0.282+0.344e-3*x 1 1.e10 s

*vec/rea a1,b1 loss.2 ! 'OC'
*vec/pl b1%a1 ! sl



nnx=9
max=[num]-11

set PLCI 4
set lwid 6
set dmod 2
vec/pl c%a ! sl

set PLCI 3
set lwid 3
set dmod 1
vec/cre c2(100) R
vec/rea a2,b2,c2 loss.2 ! 'OC'
sigma b2=b2/360
sigma c2=c2/sqrt(360)
vec/pl b2%a2 ! sl
vec/pl c2%a2 ! sl


set PLCI 7
set lwid 3
set dmod 1
sigma er=0.001*c
vec/fit a([nnx]:[max]) c([nnx]:[max]) er([nnx]:[max]) z se0 2 xpar
f/pl z 1e-1 1e11 s
vec/inp xpar(3) 360
appl comis quit
        call parout
        end
quit
f/pl z 1e-1 1e11 s

set PLCI 6
set lwid 6
set dmod 1

set PLCI 6
set lwid 6
set dmod 1
annx=$SIGMA(a([nnx]))
amax=$SIGMA(a([max]))
b1nnx=$SIGMA(c([nnx])*3.)
b2nnx=$SIGMA(c([nnx])*.3)
b1max=$SIGMA(c([max])*3.)
b2max=$SIGMA(c([max])*.3)
dline [annx] [annx] [b1nnx] [b2nnx]
dline [amax] [amax] [b1max] [b2max]

set PLCI 7
set lwid 3
set dmod 1


set PLCI 4
set lwid 7
set dmod 2
gra [num] a b ls

opt linx
opt liny
null 0 1 0 1 sa

set TXCI 3
set chhe 0.4
itx 0.65 0.20 "L#dE/dx"G#=a+bE
itx 0.60 0.15 a=$FORMAT(par1(1), f5.3) "M# GeV/mwe "N#,
itx 0.60 0.10 b=$FORMAT(par1(2)*1000, f5.3) "M# 10^-3!/mwe "N#
itx 0.10 0.90 "L#([d]E^2!/dx)"G#^1/2!=a^*!+b^*!E
itx 0.10 0.85 a^*!=$FORMAT(xpar(1), f5.3) "M# GeV/mwe^1/2! "N#,
itx 0.10 0.80 b^*!=$FORMAT(xpar(2)*1000, f5.3) "M# 10^-3!/mwe^1/2! "N#
set TXCI 1
atit 'energy "M# GeV "N#' 'total energy losses "M# GeV/mwe "N#'

exec plot1 loss
return


macro ang
exec init1
flag=2

nt/cre 10 ' ' 4 ' ' 10000 t p dt dp
nt/rea 10 scat1

nt/cre 20 ' ' 4 ' ' 10000 t2 p2 dt2 dp2
nt/rea 20 scat2

set XVAL 0.4
null -1 1 -1 1
dline -0.3 0.3 -0.7 -0.7
dline 0 0 -0.85 -0.55

dline -0.3 0.3 0 0
dline 0 0 -0.15 0.15

set MTYP 7

set ncol 58
palette 1

if [flag].eq.1 then

appl comis quit
        real function fput
        real t, p, dt, dp
        logical chain
        character*128 cfile
        common /pawchn/ chain,nchevt,ichevt
        common /pawchc/ cfile
        common /pawidn/ idnevt,obs(13), t, p, dt, dp
	dimension xx(1), yy(1)
	if ( p.gt.180 ) then
		p=360-p
	endif
	g=float(nint(p/4)+10)
	call igset('PMCI', g)
	xx(1)=dp
	yy(1)=dt
	call ipm(1, xx, yy)
        end
quit

nt/lo 10 fput

elseif [flag].eq.2 then

do i=90,1,-1

if [i].le.90 then
	n=[i]
else
	n=180-[i]
endif

set PMCI 45-[n]/2+10

min=([i]-1)
max=[i]

do j=360,10,-10

jmin=([j]-10)
jmax=[j]

nt/pl 10.dt%dp (t.ge.[min].and.t.lt.[max]).and.(p.ge.[jmin].and.p.lt.[jmax]) ! ! ! s

enddo

nt/pl 20.dt2-0.7%dp2 (t2.ge.[min].and.t2.lt.[max]) ! ! ! s

enddo

endif

atit 'deviation in [f] from direction of primary "M# [\260] "N#' 'deviation in [q] from direction of primary "M# [\260] "N#'


opt linx
opt liny
null 0 1 0 1 sa

set TXCI 4
set chhe 0.5
itx .32 .65 with magnetic field:

itx .34 .25 no magnetic field:

itx .2 .8 [q] color map:
appl comis quit
	subroutine phi
	data x/0.65/, y/0.8/

	call isfais(1)
	do i=1,180
	if ( i.le.90 ) then
		n=i
	else
		n=180-i
	endif
	fmin=float(i-1)
	fmax=float(i)
	call isfaci(n/2+10)
	call igarc(x, y, 0.0, 0.1, fmin, fmax)
	enddo

	call isfaci(55)
	call igarc(x, y, 0.0, 0.18, 89., 91.)

	call isfaci(40)
	call igarc(x, y, 0.0, 0.18, 59., 61.)

	call isfaci(1)
	call igarc(x, y, 0.121, 0.125, 62., 88.)

	end
	call phi
	end
quit

x=0.7077
y=0.9086

x1=[x]-.020
y1=[y]+.018

x2=[x]-.027
y2=[y]+.003

dline [x] [x1] [y] [y1]
dline [x] [x2] [y] [y2]

itx 0.675 0.935 [q]

exec plot1 ang
return


macro init1
set *
opt *
pic/del *
vec/del *
hi/del *
opt zfl
opt NBOX
set XMGL 3.0
set YMGL 2.4
set XVAL 0.8
set YVAL 0.4
set XLAB 2.0
set YLAB 1.4
set KSIZ .50
set GSIZ .50
set TSIZ .50
set ASIZ .50
set CSIZ .50
set PSIZ .50
set VSIZ .50
set SSIZ .50
set 2SIZ .50
igset LASI .50
igset TMSI .50
set HWID 5
set BWID 1
set PWID 5
set FWID 5
set XWID 1
set YWID 1
set CWID 5
igset LWID 5
igset TXFP -130
set CHHE .50
set cfon -40
set gfon -60
set lfon -40
set tfon -40
set vfon -40
return


macro plot1
fort/file 66 [1].ps
meta 66 -111
pict/plot *
close 66
sh rm -f [1].ps.gz
sh gzip [1].ps
sh rm -f a.ps.gz
sh ln -s [1].ps.gz a.ps.gz
return