*
*       SIMULNET.F, PROGRAM USED TO SIMULATE THE HETEROGENEOUS AGENT
*       MODELS IN MY PAPER
*
*       "SOLVING DYNAMIC MODELS WITH AGGREGATE SHOCKS AND
*        HETEROGENEOUS AGENTS"
*
*       Journal of Macroeconomic Dynamics (forthcoming)
*
*       MARCH 19, 1997
*
*       this program simulates EXP2Q30, EXP2Q500, and EXP4Q101
*
*
	implicit real*8(a-h,o-z),integer(i-n)
	parameter(npar=50,npar2=20)
	common/solut/   aold(npar),anew(npar)
	common/dat/     dat(100000,10),adat(10),sddat(10),ydat(10,10),
     $                  audat(10,10),yadat(10,10),qdat(10,10)
	common/state/   stu(1000,10),stbar(10,2),st(1000,10)
	common/rpar/    xg,blow,bhigh,bdiff,bondp,beta,gam,price(10)               
	common/ipar/    i1max,i2max,i3max,i4max,i5max,i6max,nsimpson,
     $  is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome
	common/prob/    prob(2,2),qx(24),qw(24),pagg(2,2)
	common/coef/    aggold(4,100),aggnew(4,100)
	common/hermite/ hx(1000),hw(1000)
	

	dimension bnew(100000),ww(100000),root(1000),
     $  roothx(1000),roothw(1000)
	
c
c
c       npar1 is number of parameters for individual policy function
c       npar2 is number of parameters for distribution an bond price
c

	open(36,file= 'cheb0.dat',   status='unknown')
	open(38,file= 'shock.dat',  status='unknown')
	open(31,file= 'bigcoef.b33',   status='unknown')
	open(32,file= 'bigagg.b33',  status='unknown')
	open(33,file= 'bigpar.b33',    status='unknown')
	open(331,file='bigcoef.333',   status='unknown')
	open(332,file='bigagg.333',  status='unknown')
	open(333,file='bigpar.333',    status='unknown')
	open(41,file= 'bigcoef.b44',   status='unknown')
	open(42,file= 'bigagg.b44',  status='unknown')
	open(43,file= 'bigpar.b44',    status='unknown')
	open(81,file= 'bigcoef.bi8',   status='unknown')
	open(82,file= 'bigagg.bi8',  status='unknown')
	open(83,file= 'bigpar.bi8',    status='unknown')

*
*       the "*.b33" files correspond with EXP2Q30
*       the "*.333" files correspond with EXP2Q500
*       the "*.b44" files correspond with EXP4Q101 (tight bc)
*       the "*.bi8" files correspond with EXP4Q101 (relaxed bc)
*

	read(36,*) (root(j),j=1,870)



c
c
c       prob(i,j) = probability of going from j to i
c
c


	write(78,*) 'OUTPUT NOTATION CONSISTENT WITH THE PAPER'

c
c	don't ask me why I did it like this but note that
c	when ij2 = 1 or 2 the value of ij1 doesn't matter
c

	do 66 ij1 = 1,2        
	do 66 ij2 = 1,3
	write(*,*) ij1,ij2
	idum = -84

	if(ij2.le.2) then
	open(37,file= 'pr.dat',    status='unknown')
	read(37,*) prob(1,1),prob(2,2),prob(1,2),prob(2,1),
     $  pagg(1,1),pagg(2,2),pagg(1,2),pagg(2,1)
	close(37)	
      else
	open(377,file= 'prob.dat',    status='unknown')
	read(377,*) prob(1,1),prob(2,2),prob(1,2),prob(2,1),
     $  pagg(1,1),pagg(2,2),pagg(1,2),pagg(2,1)
      close(377)	
	endif     

	if(ij2.eq.1) then
	read(33,*) i1max,i2max,i3max,i4max,i5max,i6max,
     $  is1max,is2max,is3max,is4max,is5max,is6max,
     $  blow,bhigh,bondp,beta,gam,eps1,eps2,
     $  bl4,bu4,bl5,bu5,bl6,bu6
	read(31,*) (aold(j),j=1,npar)
	read(32,*) ((aggold(i,j),j=1,npar2),i=1,4)

	endif
	if(ij2.eq.2) then
	read(333,*) i1max,i2max,i3max,i4max,i5max,i6max,
     $  is1max,is2max,is3max,is4max,is5max,is6max,
     $  blow,bhigh,bondp,beta,gam,eps1,eps2,
     $  bl4,bu4,bl5,bu5,bl6,bu6
	read(331,*) (aold(j),j=1,npar)
	read(332,*) ((aggold(i,j),j=1,npar2),i=1,4)
	endif
	if(ij1.eq.1.and.ij2.eq.3) then
	read(43,*) i1max,i2max,i3max,i4max,i5max,i6max,
     $  is1max,is2max,is3max,is4max,is5max,is6max,
     $  blow,bhigh,bondp,beta,gam,eps1,eps2,
     $  bl4,bu4,bl5,bu5,bl6,bu6
	read(41,*) (aold(j),j=1,npar)
	read(42,*) ((aggold(i,j),j=1,npar2),i=1,4)
	endif
	if(ij1.eq.2.and.ij2.eq.3) then
	close(41)
	close(42)
	close(43)
	read(83,*) i1max,i2max,i3max,i4max,i5max,i6max,
     $  is1max,is2max,is3max,is4max,is5max,is6max,
     $  blow,bhigh,bondp,beta,gam,eps1,eps2,
     $  bl4,bu4,bl5,bu5,bl6,bu6
	read(81,*) (aold(j),j=1,npar)
	read(82,*) ((aggold(i,j),j=1,npar2),i=1,4)
	endif


	imax   =   ipoint(i1max,i2max,i3max) + 18
	
*
*
*       Construct Gridpoints 
*
*
*       stbar(.,.) contains the boundaries
*
*

	ifrac1 = 0        
	ifrac2 = 0

	stu(1,2) = 0.3772
	stu(2,2) = 0.6228
	stu(1,3) = 0.991
	stu(2,3) = 1.047

	stbar(1,1) = stu(1,2) +  blow
	stbar(1,2) = stu(2,2) +  bhigh
	stbar(4,1) = bl4
	stbar(4,2) = bu4
	stbar(5,1) = bl5
	stbar(5,2) = bu5
	stbar(6,1) = bl6
	stbar(6,2) = bu6
c
c
c       note that implicitly there is a scaling for the discrete
c       valued state variables (i=2,3).
c
c       Since you never interpolate these values, you don't care
c       what this scaling is  [it is not sca(.)]

c
c
c       below we first find the appropriate roots, and then
c       find the corresponding unscaled numbers
c
c

	do 10 is1 = 1,is1max
	istart = (is1max-2)*30
	st(is1,1) = root(istart+is1)
10      continue

	do 20 is1 = 1,is1max
	bb = st(is1,1)
	stu(is1,1) = scinv(bb,1)
20      continue
		

	do 11 is2 = 1,is2max
	istart = (is2max-2)*30
	st(is2,2) = root(istart+is2)
11      continue
 
	do 12 is3 = 1,is3max
	istart = (is3max-2)*30
	st(is3,3) = root(istart+is3)
12      continue
	 
	do 13 is4 = 1,is4max
	istart = (is4max-2)*30
	st(is4,4) = root(istart+is4)
13      continue
	do 21 is4 = 1,is4max
	bb = st(is4,4)
	stu(is4,4) = scinv(bb,4)
21      continue
		
	do 14 is5 = 1,is5max
	istart = (is5max-2)*30
	st(is5,5) = root(istart+is5)
14      continue
	do 22 is5 = 1,is5max
	bb = st(is5,5)
	stu(is5,5) = scinv(bb,5)
22      continue

	do 15 is6 = 1,is6max
	istart = (is6max-2)*30
	st(is6,6) = root(istart+is6)
15      continue
	do 23 is6 = 1,is6max
	bb = st(is6,6)
	stu(is6,6) = scinv(bb,6)
23      continue
     
	

*        
*
*       start simulation across time
*
*
*
*
*       first period
*
*

*
*       define the four aggregate state variables
*       the distribution variables are for the
*       begin-of-period bondholdings relative to
*       the current aggregate endowment
*
*

	i3new =  1
	xg    = stu(i3new,3)
	s4new = stu(3,4)
	s5new = stu(3,5)
	s6new = stu(3,6)
	s4new = sca(s4new,4)
	s5new = sca(s5new,5)
	s6new = sca(s6new,6)
	s3new = st(i3new,3)
	
*
*       
*       bond price at the state variables according to the
*       numerical solution
*
*

	ieq   = 1
	bondp = aggfunc(s3new,s4new,s5new,s6new,ieq)
	
*
*       calculate the choice variables for our agent
*       for this bond price
*
*

	ww(1)  =  0.5
	inew   =  1


	s1new      = sca(ww(1),1)        
	s2new = st(inew,2)
	exp1 = pol(npar,aold,s1new,s2new,s3new,s4new,s5new,s6new)
	cons1   =   ( beta*exp1 / bondp )**(1./gam)
	bnew(1) =  ( ww(1)  - cons1 ) / bondp
	if(bnew(1).lt.blow) then
	bnew(1) =  blow
	cons1 = ww(1) - bondp*bnew(1)
	endif
	if(bnew(1).gt.bhigh) then
	bnew(1) =  bhigh
	cons1 = ww(1) - bondp*bnew(1)
	endif

910     continue


*
*
*       update what is old
*
*

	iold = inew
	cold = cons1
	s3old = s3new
	s4old = s4new
	s5old = s5new
	s6old = s6new

*
*       other period's
*
*

	itot = 100000
	
	do 448 itime = 2,itot
	
*
*       find the new draw for the aggregate growth rate
*


	r3 = ran1(idum)


	if(i3old.eq.1) then
		if(r3.lt.pagg(1,1)) then
		i3new = 1
		else
		i3new = 2
		endif
	endif
	if(i3old.eq.2) then
		if(r3.lt.pagg(2,2)) then
		i3new = 2
		else
		i3new = 1
		endif
	endif
	s3new = st(i3new,3)

*
*
*       use the numerical solution to find the new endogenous
*       aggregate state variables and the bond price
*       
*

	ieq   = 2
	s4new = aggfunc(s3new,s4new,s5new,s6new,ieq)
	ieq   = 3
	s5new = aggfunc(s3new,s4new,s5new,s6new,ieq)
	ieq   = 4
	s6new = aggfunc(s3new,s4new,s5new,s6new,ieq)

	z4    = scinv(s4new,4)
	z5    = scinv(s5new,5)
	z6    = scinv(s6new,6)

	if(ij2.le.2) then
	z4      = z4/(stu(i3new,3))
	s4new   = sca(z4,4)
	z5      = z5/(stu(i3new,3)**2.)
	s5new   = sca(z5,5)
	z6      = z6/(stu(i3new,3)**2.)
	s6new   = sca(z6,6)
	endif

	if(ij2.gt.2) then
	z4      = z4/(stu(i3new,3)**2.)
	s4new   = sca(z4,4)
	z5      = z5/(stu(i3new,3)**3.)
	s5new   = sca(z5,5)
	z6      = z6/(stu(i3new,3)**4.)
	s6new   = sca(z6,6)
	endif

	ieq   = 1
	s7    = 0.
	bondp = aggfunc(s3new,s4new,s5new,s6new,ieq)

*
*
*       find the new draws for the idiosyncratic shocks,
*       calcuate the new wealth levels
*
*


	r1 = ran1(idum)
	
	if(iold.eq.1) then
		if(r1.lt.prob(1,1)) then
		inew = 1
		else
		inew = 2
		endif
	endif

	if(iold.eq.2) then
		if(r1.lt.prob(2,2)) then
		inew = 2
		else
		inew = 1
		endif
	endif

	ww(itime)  = stu(inew,2) + bnew(itime-1) / stu(i3new,3)
	
*        
*
*       calcuate the agent's decisions
*
*

	s2new      = st(inew,2)
	s1new      = sca(ww(itime),1)        

	exp1 = pol(npar,aold,s1new,s2new,s3new,s4new,s5new,s6new)
	
	cons1       =   ( beta*exp1 / bondp )**(1./gam)
	bnew(itime) =  ( ww(itime)  - cons1 ) / bondp

	if(bnew(itime).lt.blow) then
	bnew(itime) =  blow
	cons1 = ww(itime) - bnew(itime)*bondp
	if(itime.gt.100) ifrac1 = ifrac1 + 1
	endif
	if(bnew(itime).gt.bhigh) then
	bnew(itime) =  bhigh
	cons1 = ww(itime) - bnew(itime)*bondp
	if(itime.gt.100) ifrac2 = ifrac2 + 1
	endif


	if(itime.ge.901.and.itime.le.1000) 
     $  write(79,'(i6,5f14.7)') itime,stu(inew,2),stu(i3new,3),
     $  cons1,bnew(itime),bondp

*
*
*       update what is old before you go to the next period
*
*

       

	dat(itime,1) = log(cons1)-log(cold)+log(stu(i3new,3))
	dat(itime,2) = bnew(itime)*2.
	dat(itime,3) = stu(inew,2)
	dat(itime,4) = stu(i3new,3)*2.
	dat(itime,5) = (1.-bondp)/bondp

	cold = cons1
	i3old = i3new
	s3old = s3new
	s4old = s4new
	s5old = s5new
	s6old = s6new
	iold = inew

448     continue

	write(*,'(2f14.6)') dble(ifrac1)/dble(itot-100),
     $                      dble(ifrac2)/dble(itot-100)
	write(78,'(2f14.6)') dble(ifrac1)/dble(itot-100),
     $                       dble(ifrac2)/dble(itot-100)

	do 710 j = 1,5
	adat(j)  =  ave(j,itot)
	sddat(j) = (var(j,itot))**0.5
	write(*,'(i6,2f14.6)') j,adat(j),sddat(j)
	write(78,'(i6,2f14.6)') j,adat(j),sddat(j)
710     continue

	do 740 k = 1,2
	
	do 720 l = 3,7
	audat(k,l) = aucor(k,l,itot)
	yadat(k,l) = yacor(k,l,itot)
	qdat(k,l)  = qcor(k,l,itot)
	ydat(k,l)  = ycor(k,l,itot)
720     continue
	write(*,'(5f12.6)') (audat(k,l),l=3,7)
	write(*,'(5f12.6)') (yadat(k,l),l=3,7)
	write(*,'(5f12.6)') (qdat(k,l),l=3,7)
	write(*,'(5f12.6)') (ydat(k,l),l=3,7)
	write(78,'(5f12.6)') (audat(k,l),l=3,7)
	write(78,'(5f12.6)') (yadat(k,l),l=3,7)
	write(78,'(5f12.6)') (qdat(k,l),l=3,7)
	write(78,'(5f12.6)') (ydat(k,l),l=3,7)

740     continue

66      continue


	stop
	end


	real*8 function aucor(j,l,itot)
	implicit real*8(a-h,o-z),integer(i-n)
	common/dat/     dat(100000,10),adat(10),sddat(10),ydat(10,10),
     $                  audat(10,10),yadat(10,10),qdat(10,10)
	
	aaa = ave(j,itot)
	vvv = var(j,itot)
	lag = -l+5
	cov = 0.
	do 100  i = 96,itot-5
	cov = cov + (dat(i,j)-aaa)*(dat(i-lag,j)-aaa)
100     continue
	cov = cov / dble(itot-100)
	aucor = cov/vvv

	end

	real*8 function yacor(j,l,itot)
	implicit real*8(a-h,o-z),integer(i-n)
	common/dat/     dat(100000,10),adat(10),sddat(10),ydat(10,10),
     $                  audat(10,10),yadat(10,10),qdat(10,10)
	
	aa1 = ave(j,itot)
	aa2 = ave(4,itot)
	sd1 = (var(j,itot))**0.5
	sd2 = (var(4,itot))**0.5
	lag = -l+5
	cov = 0.
	do 100  i = 96,itot-5
	cov = cov + (dat(i,j)-aa1)*(dat(i-lag,4)-aa2)
100     continue
	cov = cov / dble(itot-100)
	yacor=cov/(sd1*sd2)

	end

	real*8 function qcor(j,l,itot)
	implicit real*8(a-h,o-z),integer(i-n)
	common/dat/     dat(100000,10),adat(10),sddat(10),ydat(10,10),
     $                  audat(10,10),yadat(10,10),qdat(10,10)
	
	aa1 = ave(j,itot)
	aa2 = ave(5,itot)
	sd1 = (var(j,itot))**0.5
	sd2 = (var(5,itot))**0.5
	lag = -l+5
	cov = 0.
	do 100  i = 96,itot-5
	cov = cov + (dat(i,j)-aa1)*(dat(i-lag,5)-aa2)
100     continue
	cov = cov / dble(itot-100)
	qcor= cov/(sd1*sd2)

	end

	real*8 function ycor(j,l,itot)
	implicit real*8(a-h,o-z),integer(i-n)
	common/dat/     dat(100000,10),adat(10),sddat(10),ydat(10,10),
     $                  audat(10,10),yadat(10,10),qdat(10,10)
	
	aa1 = ave(j,itot)
	aa2 = ave(3,itot)
	sd1 = (var(j,itot))**0.5
	sd2 = (var(3,itot))**0.5
	lag = -l+5
	cov = 0.
	do 100  i = 96,itot-5
	cov = cov + (dat(i,j)-aa1)*(dat(i-lag,3)-aa2)
100     continue
	cov = cov / dble(itot-100)
	ycor= cov/(sd1*sd2)

	end

	real*8 function ave(j,itot)
	implicit real*8(a-h,o-z),integer(i-n)
	common/dat/     dat(100000,10),adat(10),sddat(10),ydat(10,10),
     $                  audat(10,10),yadat(10,10),qdat(10,10)
			
	ave = 0.
	do 50 i = 101,itot
	ave = ave + dat(i,j)
50      continue
	ave = ave/dble(itot-100)

	end

	real*8 function var(j,itot)
	implicit real*8(a-h,o-z),integer(i-n)
	common/dat/     dat(100000,10),adat(10),sddat(10),ydat(10,10),
     $                  audat(10,10),yadat(10,10),qdat(10,10)
			
	aaa = ave(j,itot)
	var = 0.
	do 50 i = 101,itot
	var = var + (dat(i,j)-aaa)**2.
50      continue
	var = var/dble(itot-100)

	end

	

	real*8 function sca(x,i)
	implicit real*8(a-h,o-z),integer(i-n)
	common/state/   stu(1000,10),stbar(10,2),st(1000,10)
	common/rpar/    xg,blow,bhigh,bdiff,bondp,beta,gam,price(10)               
	common/ipar/    i1max,i2max,i3max,i4max,i5max,i6max,nsimpson,
     $  is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome
	

	sca =  ( 2.*x -  stbar(i,2) - stbar(i,1) ) 
     $        /  (stbar(i,2)-stbar(i,1))

	end 

	real*8 function scinv(x,i)
	implicit real*8(a-h,o-z),integer(i-n)
	common/state/   stu(1000,10),stbar(10,2),st(1000,10)
	common/rpar/    xg,blow,bhigh,bdiff,bondp,beta,gam,price(10)               
	common/ipar/    i1max,i2max,i3max,i4max,i5max,i6max,nsimpson,
     $  is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome
	 

	scinv =  (    stbar(i,2) + stbar(i,1)  
     $        +    x*(stbar(i,2) - stbar(i,1))   ) / 2.

	end 

*
*
*        Define Chebyshev Polynomials
*
*

	real*8 function h(x,i)
	implicit real*8(a-h,o-z),integer(i-n)
	dimension a(200)

	a(1) = 1.
	a(2) = x

	if(i.le.1) then
	h =  a(i+1)
	goto 99
	endif

	
	if(i.gt.1) then

	do 10 j = 2,i
	a(j+1) = 2.*x*a(j) - a(j-1)
10      continue
	h =  a(i+1)

	endif

99      continue

	end


	real*8 function hh(s1,s2,s3,i1,i2,i3)
	implicit real*8(a-h,o-z),integer(i-n)
	
	hh = h(s1,i1)*h(s2,i2)*h(s3,i3)
	
	end

	real*8 function hh2(s3,s4,s5,s6,s7,i3,i4,i5,i6,i7)
	implicit real*8(a-h,o-z),integer(i-n)
	
	hh2 = h(s3,i3)*h(s4,i4)*h(s5,i5)*h(s6,i6)*h(s7,i7)
	
	end


	integer function ipoint(i1,i2,i3)
	implicit real*8(a-h,o-z),integer(i-n)
	common/ipar/    i1max,i2max,i3max,i4max,i5max,i6max,nsimpson,
     $  is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome
	 

	ipoint = 1 + i2 
     $  + i3*(i2max+1)  
     $  + i1*(i3max+1)*(i2max+1) 

c
c       i* goes from 0 to i*max
c
c       with this convention increasing i1max adds parameters at the
c       end of the parameter vector
c
	
	end




	real*8 function aggfunc(ss3,ss4,ss5,ss6,ireg)                 
	Implicit real*8(a-h,o-z),integer(i-n)
	parameter(npar2=20)
	common/state/   stu(1000,10),stbar(10,2),st(1000,10)
	common/rpar/    xg,blow,bhigh,bdiff,bondp,beta,gam,price(10)               
	common/ipar/    i1max,i2max,i3max,i4max,i5max,i6max,nsimpson,
     $  is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome
	common/prob/    prob(2,2),qx(24),qw(24),pagg(2,2)
	common/coef/    aggold(4,100),aggnew(4,100)
	common/hermite/ hx(1000),hw(1000)

	aggfunc = 0.
	
	do 44 i = 1,npar2
	if(i.eq.1)  term = 1
	if(i.eq.2)  term = ss3
	if(i.eq.3)  term = ss4
	if(i.eq.4)  term = ss5
	if(i.eq.5)  term = ss6
	if(i.eq.6)  term = ss4*ss3
	if(i.eq.7)  term = ss5*ss3
	if(i.eq.8)  term = ss6*ss3
	if(i.eq.9)  term = h(ss4,2) 
	if(i.eq.10) term = h(ss5,2)
	if(i.eq.11) term = h(ss6,2)
	if(i.eq.12) term = ss4*ss5
	if(i.eq.13) term = ss4*ss6
	if(i.eq.14) term = ss5*ss6
	if(i.eq.15) term = h(ss4,2)*ss3 
	if(i.eq.16) term = h(ss5,2)*ss3 
	if(i.eq.17) term = h(ss6,2)*ss3 
	if(i.eq.18) term = ss4*ss5*ss3 
	if(i.eq.19) term = ss4*ss6*ss3 
	if(i.eq.20) term = ss5*ss6*ss3 
c        write(*,*) i,term,aggold(ireg,i)
	aggfunc = aggfunc + aggold(ireg,i)*term
44      continue

	end

	real*8 function pol(npar,coef,s1,s2,s3,s4,s5,s6)
	implicit real*8(a-h,o-z),integer(i-n)
	common/state/ stu(1000,10),stbar(10,2),st(1000,10)
	common/rpar/    xg,blow,bhigh,bdiff,bondp,beta,gam,price(10)               
	common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson,
     $  is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome
     
	dimension coef(npar)
	
	pol = 0.
	do 20 j1 = 1,i1max+1
	do 20 j2 = 1,i2max+1
	do 20 j3 = 1,i3max+1
	i1 = j1 - 1
	i2 = j2 - 1
	i3 = j3 - 1
	jj   =  ipoint(i1,i2,i3)
	pol =  pol + coef(jj)*hh(s1,s2,s3,i1,i2,i3)
20      continue
	jjmax = ipoint(i1max,i2max,i3max)
	pol = pol +  coef(jjmax+1)*s4
	pol = pol +  coef(jjmax+2)*s5
	pol = pol +  coef(jjmax+3)*s6
	pol = pol +  coef(jjmax+4 )*h(s4,2)
	pol = pol +  coef(jjmax+5 )*h(s5,2)
	pol = pol +  coef(jjmax+6 )*h(s6,2)
	pol = pol +  coef(jjmax+7 )*s4*s5
	pol = pol +  coef(jjmax+8 )*s4*s6
	pol = pol +  coef(jjmax+9 )*s5*s6
	pol = pol +  coef(jjmax+10)*s4*s3
	pol = pol +  coef(jjmax+11)*s5*s3
	pol = pol +  coef(jjmax+12)*s6*s3
	pol = pol +  coef(jjmax+13)*h(s4,2)*s3
	pol = pol +  coef(jjmax+14)*h(s5,2)*s3
	pol = pol +  coef(jjmax+15)*h(s6,2)*s3
	pol = pol +  coef(jjmax+16)*s4*s5*s3
	pol = pol +  coef(jjmax+17)*s4*s6*s3
	pol = pol +  coef(jjmax+18)*s5*s6*s3
	end

	REAL*8 FUNCTION RAN1(IDUM)
c
c       returns a uniform random deviae between 0.0 and 1.0.
c       set IDUM to any negative value to initialize or reinitialize
c       the sequence

	IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N)

	DIMENSION R(97)
	PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=1./M1)
	PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=1./M2)
	PARAMETER (M3=243000,IA3=4561,IC3=51349)
	save iff,ix1,ix2,ix3,r
	DATA IFF /0/
	IF(IDUM.LT.0.OR.IFF.EQ.0) THEN
		IFF = 1
		IX1 = MOD(IC1-IDUM,M1)
		IX1 = MOD(IA1*IX1+IC1,M1)
		IX2 = MOD(IX1,M2)
		IX1 = MOD(IA1*IX1+IC1,M1)
		IX3 = MOD(IX1,M3)
		DO 11 J=1,97
			IX1=MOD(IA1*IX1+IC1,M1)
			IX2=MOD(IA2*IX2+IC2,M2)
			R(J)=(DBLE(IX1)*DBLE(IX2)*RM2)*RM1
11              CONTINUE

		IDUM = 1

	ENDIF

	IX1=MOD(IA1*IX1+IC1,M1)
	IX2=MOD(IA2*IX2+IC2,M2)
	IX3=MOD(IA3*IX3+IC3,M3)
	J=1+(97*IX3)/M3
	IF(J.GT.97.OR.J.LT.1)PAUSE
	RAN1=R(J)
	R(J)=(DBLE(IX1)+DBLE(IX2)*RM2)*RM1
	RETURN
	END


