module iomod
	use likmod
	
	implicit none
	
	contains
	
	subroutine plotspecdens(th)
		real				:: th(nth)
		integer, parameter	:: npoints=200
		real, parameter		:: lammax=hub+10
		integer				:: i
		real, parameter		:: lamgrid(npoints)=[(lammax*(i-1.0)/npoints,i=1,npoints)]
		real, allocatable, save	:: mp(:,:)
		integer, save		:: n=0
				
		if(n==0) allocate(mp(0,0))
		n=min(n+1,7)
		mp=reshape([getldens(th),mp(:,1:n-1)],[npoints,n])
		call storeML(mp,"y")
		call storeML(lamgrid,"x")
		call execinML("plot(x,y)")
		
	contains
	
		function getldens(th) result(val)
			use polymod
			real		:: th(nth),val(npoints)
			real		:: apoly(p+1), bpoly(p), aapoly(2*p+1), bbpoly(2*p-1)
			integer		:: i
			apoly=realpolyfromroots(-getcfromh(th(1:p)))
			aapoly=polymult(apoly,polyminus(apoly))
			if(p==1) then
				bpoly=1
				bbpoly=1
			else
				bpoly=realpolyfromroots(-getcfromh(th(p+1:)))
				bbpoly=polymult(bpoly,polyminus(bpoly))
			endif
			do i=1,npoints
				val(i)=polyevalcmplx(bbpoly,cmplx(0,lamgrid(i)))/polyevalcmplx(aapoly,cmplx(0,lamgrid(i)))
			enddo
!			val=log(val*(para%c0**2+lamgrid**2))
			val=log(val)
		end function
	end subroutine
end module