c   pgdists.f:

c   To compile:

c   f77 -O -o pgdistsx pgdists.f dists.o -L/usr/local/lib\
c              -L/usr/openwin/lib -lpgplot -lX11

c   To run:

c   #!/bin/tcsh
c   setenv PGPLOT_DIR /usr/local/lib/pgplot
c   setenv PGPLOT_DEV /xwin
c   exec ile -/usr/local/Splus3.4/.ilerc pgdistsx

	double precision zpdf,tpdf,chipdf,fpdf
	double precision dx,dcdf,dpdf,dq,du1,du2
	character*40 gname,xlab,ylab,mlab,dflab

	dimension x(1000),y(1000)
	integer pgbegin

	write(*,1)
 1	format(' Enter nx,ny: '$)
	read(*,*) nx,ny

	write(*,2)
 2	format(' Enter graphics device: ',/,
     1         '''/xwin'' or ''fname/ps'' or ''fname/vps'': '$)
	read(*,*) gname

	if(pgbegin(0,gname,nx,ny).ne.1) stop

	u1=.01
	u2=.99
	npts=51


  5   	write(*,10)
 10   	format(' 1: Z',/,
     1         ' 2: t',/,
     1         ' 3: chi-square',/,
     1         ' 4: F',/,
     1         ' 5: Change u1,u2',/,
     1         ' 6: STOP')

      	read(*,*) iopt
      	if(iopt.lt.1.or.iopt.gt.6) go to 5
      
  	if(iopt.eq.6) go to 99

	if(iopt.eq.5) then
	
  16		write(*,20)
  20		format(' Enter u1, u2: ',$)
		read(*,*) u1,u2
		if(u1.lt.0.0.or.u1.gt.1.0.or.u1.gt.u2) go to 16
		go to 5
	endif

 25   	write(*,30)
 30   	format(' 1: pdf',/,
     1         ' 2: cdf',/,
     1         ' 3: quantile',/,
     1         ' 4: BACK TO MAIN MENU')

      	read(*,*) iopt1
	if(iopt1.lt.0.or.iopt1.gt.4) go to 25

	if(iopt1.eq.4) go to 5

	if(iopt.eq.2.or.iopt.eq.3) then
  35		write(*,40)
  40		format(' Enter df: ',$)
		read(*,*) nd1
		if(nd1.lt.1) go to 35
	endif

	if(iopt.eq.4) then
  45		write(*,50)
  50		format(' Enter df1, df2: ',$)
		read(*,*) nd1,nd2
		if(nd1.lt.0.or.nd2.lt.0) go to 45
	endif

	umin=u1
	umax=u2

	if(iopt1.eq.1.or.iopt1.eq.2) then

		if(iopt.eq.1) then
			call zqnt(dble(u1),du1,ier)
			call zqnt(dble(u2),du2,ier)
		endif
		if(iopt.eq.2) then
			call tqnt(dble(u1),nd1,du1,ier)
			call tqnt(dble(u2),nd1,du2,ier)
		endif
		if(iopt.eq.3) then
			call chiqnt(dble(u1),nd1,du1,ier)
			call chiqnt(dble(u2),nd1,du2,ier)
		endif
		if(iopt.eq.4) then
			call fqnt(dble(u1),nd1,nd2,du1,ier)
			call fqnt(dble(u2),nd1,nd2,du2,ier)
		endif

		umin=du1
		umax=du2
	endif


	del=(umax-umin)/(npts-1)
	do 100 i=1,npts
 100	x(i)=umin+(i-1)*del		

	do 120 i=1,npts

		dx=dble(x(i))
	
		if(iopt1.eq.1) then
			if(iopt.eq.1) y(i)=zpdf(dx)
			if(iopt.eq.2) y(i)=tpdf(dx,nd1)
			if(iopt.eq.3) y(i)=chipdf(dx,nd1)
			if(iopt.eq.4) y(i)=fpdf(dx,nd1,nd2)
		endif

		if(iopt1.eq.2) then
			if(iopt.eq.1) call zcdf(dx,dcdf,dpdf)
			if(iopt.eq.2) call tcdf(dx,nd1,dcdf,dpdf)
			if(iopt.eq.3) call chicdf(dx,nd1,dcdf,dpdf)
			if(iopt.eq.4) call fcdf(dx,nd1,nd2,dcdf,dpdf)
			y(i)=dcdf
		endif

		if(iopt1.eq.3) then
			if(iopt.eq.1) call zqnt(dx,dq,ier)
			if(iopt.eq.2) call tqnt(dx,nd1,dq,ier)
			if(iopt.eq.3) call chiqnt(dx,nd1,dq,ier)
			if(iopt.eq.4) call fqnt(dx,nd1,nd2,dq,ier)
			y(i)=dq
		endif		

 120	continue

	ymin=y(1)
	ymax=y(1)
	do 125 i=2,npts
        if(y(i).lt.ymin) ymin=y(i)
 125	if(y(i).gt.ymax) ymax=y(i)

	call pgrnge(umin,umax,umin1,umax1)
	call pgrnge(ymin,ymax,ymin1,ymax1)

	call pgenv(umin1,umax1,ymin1,ymax1,0,1)
	call pgline(npts,x,y)
	call pgpt(npts,x,y,-3)

	xlab='u'
	if(iopt1.le.2) xlab='x'
	ylab='f(x)'
	if(iopt1.eq.2) ylab='F(x)'
	if(iopt1.eq.3) ylab='Q(u)'

	mlab='Standard Normal Distribution'
	if(iopt.eq.2) mlab='Student t Distribution'
	if(iopt.eq.3) mlab='Chi Square Distribution'
	if(iopt.eq.4) mlab='F Distribution'

	call pglab(xlab,ylab,mlab)

	dflab=' '
	if((iopt.eq.2).or.(iopt.eq.3)) write(dflab,130) nd1
 130	format('df = ',i3)
	if(iopt.eq.4) write(dflab,140) nd1,nd2
 140	format('df = ',i3,1x,i3)

	call pgptxt(umin,ymax,0.0,0.0,dflab)

	go to 5

  99	continue
	call pgend

	stop
	end