module priormod
	use basicmodule
	use polymod
	implicit none
	
	
	contains
	
	subroutine sethhprior0(hhdraws)			! store (estimated) prior density on half-life bins so that we can divide by it later
		real		:: hhdraws(:)
		integer	:: n,l,i
		n=size(hhdraws)
		hhprior0(:,p)=0
		do l=1,n
			if(hhdraws(l)<=minval(hhgrid) .or. hhdraws(l)>maxval(hhgrid)) cycle
			i=count(hhdraws(l)>hhgrid)
			hhprior0(i,p)=hhprior0(i,p)+1
		enddo
		hhprior0(:,p)=nhhbins*hhprior0(:,p)/n
		hhpriorset=.true.
		if(p==pmax) then
			print *,"hhprior0"
			call mdisp(hhprior0(:,1:p))
		endif
	end subroutine
	
	function defaulthhprior(hh) result(val)		! prior before adjustment; non-uniform to ensure enough mass for all half-lives to be able to compute adjusted prior
		real		:: hh,val
		if(hh>hhub) then
			val=0
		else
			val=max(hh,hhlb)**2
			if(priorkappa>500 .and. p<3) val=(hh+10)**(5-p)
		endif
	end function
	
	function gethhprior(hh) result(val)		! adjusted prior
		real		:: hh, val
		integer	:: i
		real, parameter	:: midgrid(nhhbins)=[(0.5*(hhgrid(i-1)+hhgrid(i)),i=1,nhhbins)]

		if(.not. hhpriorset) then
			val=defaulthhprior(hh)
			return
		endif
		if(hh>hhub) then
			val=defaulthhprior(hh)
			return
		endif
		i=count(hh>hhgrid)
		val=defaulthhprior(hh)/hhprior0(i,p)	! divided so that resulting prior is uniform on half-lives
	end function

	function getsmoothprior(para) result(val)
		type(tpara)	:: para
		real		:: val
		complex		:: glist(p-1),clist(p)
		integer		:: i,j	
		clist=getcfromh(para%hAR(1:p))
		if(p>1) glist=getcfromh(para%hMA(1:p-1))
		val=0
		do i=1,p
			do j=1,p
				val=val+(clist(i)+clist(j))**-3
			enddo
		enddo
		do i=1,p
			do j=1,p-1
				val=val-2*(clist(i)+glist(j))**-3
			enddo
		enddo
		do i=1,p-1
			do j=1,p-1
				val=val+(glist(i)+glist(j))**-3
			enddo
		enddo
		val=-priorkappa*val
	end function		

	
	function rootdist(para) result(val)
		type(tpara)	:: para
		real		:: val
		complex	:: dlist(p-1),clist(p)
		integer		:: i,j
		if(p==1) then
			val=1E100
			return
		endif
		clist=getcfromh(para%hAR(1:p))
		val=0
		do i=1,p
			do j=1,i
				if (abs((1.0-clist(i)/Tstar)*(1.0-clist(j)/Tstar))>0.999) return
			enddo
		enddo
		val=1E100
	end function	

end module