c		c:\drw2\simple\2lenda.f
c
c
c		program to calculate:
c
c		(a) stuff for fixed H when you have two lenders
c
c		This program:
c
c		- uses simple two point distribution {0,hu}
c		- the fraction that receives hu, gamh, depends on H
c
c
c
c		December 2 1998



		use msimsl
		implicit real*8(a-h,o-z),integer(i-n)
		common/par/xvalue,alpha,prob,depr,bbbb,xi,rox,gam,eta,
     $	hstst,vcost,ststg,ststw,ststge,ubar
		common/iiii/iend,ifunc


		external fgam


		open(99,file='for099.dat',
     $	status='unknown')


		iprod = 1

		alpha = 0.9
		bbbb  = 0.96
		xi    = 0.25

		xvalue = 3.
		zero   = 0.0
		vcost  = 0.342


		ifunc   = 3
		gam     = 0.400925
		eta     = 0.2
		
		do 448 ih = 1,180

		hstst = 80.-0.5*dble(ih-1)

		ststgam  = fgam(hstst)
		ubar1    = 1.-ststgam
		ubar2    = (1.-ststgam)**2.
		ststhu   = hstst/ststgam
		ststge1  = (1.-ubar1)*bbbb*xvalue
		ststge2  = (1.-ubar2)*bbbb*xvalue

		expv1    = ststgam*ststhu**alpha

		if(iprod.eq.1) then
		expv2    = ststgam**2.*(2.*ststhu)**alpha +
     $    2.*(1-ststgam)*ststgam*ststhu**alpha
		endif
		
		if(iprod.eq.2) then
		expv2    = ststgam**2.*2.*(ststhu)**alpha +
     $    2.*(1-ststgam)*ststgam*ststhu**alpha
		endif

		ratuv    = (vcost/(ststge1*xi))**2.

		probe    = xi*ratuv**0.5
		probl    = xi/(ratuv**0.5)

		ratn1v   = probe/(ubar1+2.*probe)
		ratn2v   = probe*2.*ratn1v/ubar2
		rat1v    = ratuv+(1.-ubar1)*ratn1v+(1.-ubar2)*ratn2v
		
		ststv    = 1./rat1v
		ststn1   = ratn1v*ststv
		ststn2	 = ratn2v*ststv
		ststu    = 1.-ststn1-ststn2+ubar1*ststn1+ubar2*ststn2

		cg2		= bbbb*expv2/(1.-bbbb*(1.-ubar2))
		sg2		= bbbb*ubar2/(1.-bbbb*(1.-ubar2))

		denom   = 1.-bbbb*(1.-ubar1)+bbbb*(1.-ubar1)*probe

		cg1		= bbbb*expv1+bbbb*(1.-ubar1)*probe*cg2
		cg1		= cg1/denom
		sg1		= bbbb*(1.-ubar1)*probe*(sg2-bbbb)+bbbb*ubar1
		sg1		= sg1/denom

		denom   = 1.-probl*(1.-ratn1v)*sg1-
     $			  bbbb*(probl*ratn1v+1.-probl)

		ststw	= probl*(1.-ratn1v)*(cg1-ststge1)/denom

		ststg1	= cg1 + sg1*ststw
		ststg2  = cg2 + sg2*ststw


c
c
c		check some stuff
c
c


		frag1 = ststhu**alpha+ststg1-xvalue-ststw
		frag2 = ststhu**alpha+ststg2-xvalue-ststw
		xliq1 = ststhu**alpha+ststhu+ststge1-xvalue
		xliq2 = ststhu**alpha+ststhu+ststge2-xvalue

		write(*,'(8f10.4)') ststg2,ststg1,ststw,probe,probl,
     $    ratn1v,ststge1
		tttt = ststg2-ststg1-bbbb*ststw
		write(*,'(5f10.4)') frag1,frag2,xliq1,xliq2,tttt
		if(probe*(ststg2-ststg1-bbbb*ststw).lt.0.) stop
		if(frag1.lt.0.) stop


c
c		calculate return
c

		ststm1 = xvalue-ststge1
		ststm2 = xvalue-ststge2

		if(iprod.eq.1) then
		ret = ststn1*(1.-ubar1)*(ststhu**alpha-ststm1)+
     $		  0.5*ststn2*(ststgam**2.*((2.*ststhu)**alpha-ststm2)+
     $		          ststgam*(1.-ststgam)*(ststhu**alpha-ststm2))
		ret = ret/hstst
		endif

		if(iprod.eq.2) then
		ret = ststn1*(1.-ubar1)*(ststhu**alpha-ststm1)+
     $		  0.5*ststn2*(ststgam**2.*(2.*ststhu**alpha-ststm2)+
     $		          ststgam*(1.-ststgam)*(ststhu**alpha-ststm2))
		ret = ret/hstst
		endif


	write(*,'(8f9.5)') hstst,
     $ret,probe,probl,ststn1,ststn2,ststv,ststu
	write(99,'(12f9.5)') hstst,
     $ret,probe,probl,ststn1,ststn2,ststv,ubar1,ubar2,
     $ststn1*(1.-ubar1)*ststhu**alpha
     $+0.5*ststn2*(1.-ubar2)*(2.*ststhu)**alpha
	pause
448	continue
	stop
	end


		real*8 function fgam(hhhh)
		implicit real*8(a-h,o-z),integer(i-n)
		common/par/xvalue,alpha,prob,depr,bbbb,xi,rox,gam,eta,
     $          hstst,vcost,ststg,ststw,ststge,ubar
		common/iiii/iend,ifunc


		if(ifunc.eq.1) fgam = (gam*hhhh)**eta/(1.+(gam*hhhh)**eta)
		if(ifunc.eq.2) fgam = 1./(1.+exp(-gam*(hhhh-eta)))
		if(ifunc.eq.3) fgam = dmin1(gam*hhhh**eta,dble(0.999))
		
		return
		end