c	SUBROUTINE inv_entr
C
c	procedura per invertire un sistema di equazioni
c	integrali per dedurre la DEM 
c	con una sequenza iterativa basata sulla minimizzazione
c	dell'entropia
c	 si aggiungono come parametri anche i mesh points
c
c	gennaio 1992
c	si correggono i valori di demm secondo il rapporto
c	
C	marzo 1994
C	aggiornato per IDL
c
	subroutine inv_entr(all,yoss,nn,err,verr)
	dimension ytt(100),demmt1(100),sert(100),per(100)
	dimension cf1(100),all(100),yoss(100),rip(100),bestfi(100)
	dimension yint(100),yteor(100),yincr(100),deriv(100)
	dimension xint(100),verr(100),err(100),cf0(100),y2(100),hcf0(100)
	dimension passoc(100),demm1(100),yt(100)
	dimension xxint(1000)
	common/demmesh/ xt0(100),hhcf0(100)
	common/cfff/nnodi,hcf(100),xt(100)
	common/line1/line,kkt
	common/tt/ain0,ddt,mmt,att(100)
	common/ion/ela(100),ion(100),tmax(100),tttmax(100),demm(100)
	common/demmax/demmax
	common chi2,rms1,rms
	common/dev1/dev1
	common/ngauss/ngauss
	common/nsert/nsert(100)
	common/rip/rip(100)
	character*50,tttest,xttest
	character ela*5
	character ion*5
	character*1 yesno
	character*80 title
110	format(a1)
	ngauss=4
	write(*,*)'do you wish that the mesh points are parameters?[y/n]'
	call flush(6)
	read(*,'(a)')yesno
	i2=1
	if(yesno.eq.'y')i2=2
	chi0=1.e30
	amt=mmt
	atfin=ain0+kkt*ddt
	atin=ain0
c
c	si inizializza il vettore incrementi
c
	bestchi=1.e30
	do n=1,nnodi
	cf0(n)=hcf(n)
	if(i2.eq.2)cf0(n+nnodi)=xt(n)
	hcf0(n)=hcf(n)
	xt0(n)=xt(n)
	hhcf0(n)=hcf(n)
	if(i2.eq.2)hcf0(n+nnodi)=xt(n)
	end do
	write(9,*)(xt(n),n=1,nnodi)
	write(9,*)(hcf(n),n=1,nnodi)
	write(4,*)' mesh points:'
	write(4,*)(xt(n),n=1,nnodi)
	write(4,*)(hcf(n),n=1,nnodi)
	do iline=1,nn
	line=iline
	xx=all(iline)
	call ordina(xt,hcf,nnodi)
	yint(iline)= fct()
	sert(iline)=nsert(iline)
	write(9,11111)sert(iline),ela(iline),ion(iline),all(iline),
     +	yoss(iline),
     +	yint(iline),verr(iline)
	end do
C
	do k=1,i2*nnodi
	sumff=sumff+cf0(k)**2
	end do
	write(*,*)'do you want normalize the value in the mesh points?'
	call flush(6)
	read(*,'(a)')yesno
	if(yesno.ne.'y') go to 1002
c
c	si chiamma ttmax per sapere il max del prodotto g(T)*fi(T)
c
	do line=1,nn
	tttmax(line)=ttmax()
	yt(line)=tttmax(line)
	demm(line)=demmax
	rip(line)=yint(line)/yoss(line)
	end do
c
c	si trovano i valori di dem in tmax con la spline
c
	do line=1,nn
	yt(line)=tmax(line)
	end do
	do line=1,nn
	xx=yt(line)
	call spline(xt,hcf,nnodi,2e30,2.e30,y2)
	call splint(xt,hcf,y2,nnodi,xx,fy)
	demm(line)=fy
	rip(line)=yint(line)/yoss(line)
	end do
c
c	gennaio 1992
C
c	si correggono i valori di demm secondo il rapporto
c	iteor/ioss
c
	do line=1,nn
	dddd=alog10(yoss(line))-alog10(yint(line))
	if(dddd.lt.0)demm1(line)=demm(line)-abs(dddd)
	if(dddd.gt.0)demm1(line)=abs(dddd)+demm(line)
	end do
c
c	si chiama la spline per conoscere i valori nei mesh points
c
	call ordina(yt,demm1,nn)
	i=0
	ki=0
11	i=i+1
	k=0
	if(i.gt.nn)go to 12
	d1=0
	do j=i,nn
	if(yt(i).eq.yt(j))then
	k=k+1
	d1=d1+10**demm1(j)
	yyyy=yt(j)
	end if
	end do
	i=i+k-1
	ki=ki+1
	ytt(ki)=yt(i)
	demmt1(ki)=demm1(i)
	if(k.ne.0)then
	ytt(ki)=yyyy
	demmt1(ki)=alog10(d1/k)
	end if
	go to 11
12	continue
	yt(1)=xt(1)
	demm1(1)=hcf(1)
	do i=2,ki+1
	demm1(i)=demmt1(i-1)
	yt(i)=ytt(i-1)
	end do
	yt(ki+2)=xt(nnodi)
	demm1(ki+2)=hcf(nnodi)
	write(9,*)(yt(i),i=1,ki+2)
	write(9,*)(demm1(i),i=1,ki+2)

	nnnnn=ki+2

	do k=2,nnodi-1
	xx=xt(k)

	call spline(yt,demm1,nnnnn,2e30,2.e30,y2)
	call splint(yt,demm1,y2,nnnnn,xx,fy)
	cf0(k)=fy
	hcf(k)=cf0(k)
	hcf0(k)=cf0(k)
	end do
c
	write(9,*)(xt(n),n=1,nnodi)
	write(9,*)(hcf(n),n=1,nnodi)
	write(4,*)' mesh points:'
	write(4,*)(xt(n),n=1,nnodi)
	write(4,*)(hcf(n),n=1,nnodi)


	npassi=0
1002	alamd=0.
	nnpass=1

	chi2=1e30
	chi0=1.e30
1001	continue
	if(npassi.gt.0)then
	write(*,*)'enter lambda and #steps'
	call flush(6)
	read(*,*)alamd,nnpassi
	ipassi=0
	endif

	do k=1,nnodi
	cf1(k)=cf0(k)
	if(i2.eq.2)cf1(k+nnodi)=cf0(k+nnodi)
	if(i2.eq.2)hcf0(k+nnodi)=cf0(k+nnodi)
	end do
200	npassi=npassi+1
	ipassi=ipassi+1
	if(npassi.gt.1)then
	if(chi2.lt.bestchi)then
	do k=1,nnodi
	if(i2.eq.2)bestfi(k+nnodi)=cf0(k+nnodi)
	bestfi(k)=cf0(k)
	end do
	bestchi=chi2
	end if
	end if

	sumlam=0.
	sumcf=0.
	summax=0.

		do k=1,i2*nnodi
	passoc(k)=cf0(k)*.01
	cf1(k)=cf0(k)+passoc(k)
	if(k.eq.(i2*nnodi))cf1(k)=cf0(k)
	if(k.eq.(1+nnodi))cf1(k)=cf0(k)
	sum=0.

		do iline=1,nn
	line=iline
	xx=all(iline)
	tttest='   '
			do n=1,nnodi
	if(i2.eq.2)xt(n)=cf0(n+nnodi)
			if(xt(n).gt.xt(nnodi))then
	ttest='xt(n) gt xt(nnodi): cambia lambda'
	go to 111
			endif
	hcf(n)=cf0(n)
			end do
	call ordina(xt,hcf,nnodi)
	yteor(iline)= fct()
			do n=1,nnodi
	hcf(n)=cf1(n)
	if(i2.eq.2)xt(n)=cf1(n+nnodi)
			if(xt(n).gt.xt(nnodi))then
	ttest='xt(n) gt xt(nnodi): cambia lambda'
	go to 111
			endif
			end do
	call ordina(xt,hcf,nnodi)
	yincr(iline)= fct()
	deriv(iline)=(yincr(iline)-yteor(iline))/passoc(k)
	sigm2=err(iline)**2
	ssss=(yteor(iline)-yoss(iline))*deriv(iline)/sigm2
	sum=sum+ssss
		end do
c	peso=(cf0(k)**2/sumff)**.1
	tt=cf0(k)-alamd*sum
c	tt=cf0(k)-alamd*sum*peso
	if(abs(sum).gt.summax)summax=abs(sum)
	sumlam=sumlam+abs(sum)
	sumcf=sumcf+cf0(k)
	sumff=sumff+cf0(k)**2
	xttest='    '
	if(abs(tt).gt.30)then
	xttest=' tt >30; cambia lambda'
	do iline=1,nn
	end do
	do ii=1,i2*nnodi
	cf0(ii)=hcf0(ii)
	end do
	go to 111
	end if
	hcf0(k)=hcf(k)
	cf0(k)=(tt)
	cf1(k)=(tt)
	hcf(k)=tt
		end do

	do kk=1,nnodi
	if(i2.eq.2)xt(kk)=hcf(kk+nnodi)
	end do

c	calcolo chi square
c
2000	chi2=0.
	chi22=0.
	rms1=0
	rms2=0
	rms=0.
	do iline=1,nn
	line=iline
	xx=all(iline)
	call ordina(xt,hcf,nnodi)
	yint(iline)= fct()
	chi2=chi2+((yint(iline)-yoss(iline))/err(iline))**2
	rms=rms+((yint(iline)-yoss(iline))/yoss(iline))**2
	end do
	rms=sqrt(rms/(nn-1))
	ratio=sumcf/sumlam
	ratio1=sumcf/summax
	alamd1=(0.0001*ratio/nn)
	alamd2=(0.0001*ratio1/nn)
c	write(*,*) ' lambda =',alamd,' lambda1 = ',alamd1,' ratio=',ratio
c	write(*,*)'npassi=',npassi,'chi=',chi2/nn,'S.D.=',rms
c	write(*,*)(hcf(k),k=1,nnodi)
c	write(*,*)(xt(k),k=1,nnodi)

C
C	Test on the CHI2
C
               if(chi2.gt.chi0) then
c        write(*,*)' chi0=',chi0/nn,' chi1=',chi2/nn,'   S.D. = ',
c     +  rms
        write(4,*)' chi0=',chi0/nn,' chi1=',chi2/nn,'   S.D. = ',
     +  rms
        write(4,*)' lambda = ',alamd
        alamd=alamd/2.
	chi0=chi2
c
c	repeat the procedure with lambda=lambda/2
c
	go to 200
                end if
        if(chi0.lt.bestchi)then
        do k=1,i2*nnodi
        bestfi(k)=hcf0(k)
        end do
        write(*,*)bestchi/nn,chi0/nn
        write(*,*)(bestfi(k),k=1,nnodi)
        bestchi=chi0
        end if
        chi0=chi2
	if(npassi.eq.1)then
c
c	go to read new lambda
c
	go to  1001
	endif


99991	continue

c
c	repeat the procedure nnpassi times
c
	if(ipassi.lt.nnpassi)go to 200

c
	if(ipassi.ge.nnpass)then
	do iline=1,nn
	per(iline)=(yoss(iline)-yint(iline))/verr(iline)
	rip(iline)=yint(iline)/yoss(iline)
11111	format(i5,2a5,f10.3,3(1pe10.3),3(0pf10.3))
	sert(iline)=nsert(iline)
c	write(*,11111)nsert(iline),ela(iline),ion(iline),all(iline),
c     +	yoss(iline),yint(iline),verr(iline),per(iline),rip(iline),
c     &	tmax(iline)
	write(4,11111)nsert(iline),ela(iline),ion(iline),all(iline),
     +	yoss(iline),yint(iline),verr(iline),per(iline),rip(iline),
     &	tmax(iline)
	end do
c	write(*,*)'chi2 = ',chi2/nn,' S.D. = ',rms,' lambda = ',alamd,
c     &	'bestchi= ',bestchi/nn
	write(*,*)chi2/nn,rms,alamd,bestchi/nn
	call flush(6)
	write(4,*)' chi2 = ',chi2/nn
	write(4,*)' S.D. = ',rms
	write(4,*)' lambda = ',alamd
c
c	test per RMS<.01
c
	if(rms.lt.0.01)go to 100


c
	write(*,*)'do you want change lambda? y/n'
	call flush(6)
	read(*,'(a)')yesno
	if(yesno.eq.'y')then
c	write(*,*)'enter lambda'
c	read(*,*)alamd
	ipassi=0
c
c	go to read new lambda
c
	go to 1001
	end if
	end if

c
c	si plottano i risultati e si scrivono su disco
c
c
100	do iline=1,nn
	per(iline)=(yoss(iline)-yint(iline))/verr(iline)
	rip(iline)=yint(iline)/yoss(iline)
	sert(iline)=nsert(iline)
	write(9,11111)nsert(iline),ela(iline),ion(iline),all(iline),
     +	yoss(iline),
     +	yint(iline),verr(iline),per(iline),rip(iline),tmax(iline)
	write(4,11111)nsert(iline),ela(iline),ion(iline),all(iline),
     +	yoss(iline),
     +	yint(iline),verr(iline),per(iline),rip(iline),tmax(iline)
	end do
c
c	si prepara output per IDL

	write(*,*)(nsert(iline),iline=1,nn)
	write(*,'(a)')(ela(iline),iline=1,nn)
	write(*,'(a)')(ion(iline),iline=1,nn)
	write(*,*)(all(iline),iline=1,nn)
	write(*,*)(yoss(iline),iline=1,nn)
	write(*,*)(yint(iline),iline=1,nn)
	write(*,*)(verr(iline),iline=1,nn)
	write(*,*)(per(iline),iline=1,nn)
	write(*,*)(rip(iline),iline=1,nn)
	write(*,*)(tmax(iline),iline=1,nn)
	call flush(6)
c
	do k=1,i2*nnodi
	cf0(k)=hcf0(k)    	
	end do


111	continue
c
c	si scrive su un file la soluzione migliore
c
	write(18,*)nnodi
	do i=1,nnodi
	write(18,*)xt(i),bestfi(i)
	end do
	close (18)

	write(*,888)alamd,tttest,xttest
888	format(e10.4,2a)
	call flush(6)
	write(*,*)'do you wish to start from the best solution found
     +	till now?? [y]/n'
	call flush(6)
	read(*,'(a)')yesno
	nnpass=1
	if(yesno.eq.'y')then
c	write(*,*)bestchi/nn
c	write(*,*)(bestfi(i),i=1,nnodi)
	do k=1,i2*nnodi
	cf0(k)=bestfi(k)
	end do
	end if
c	write(*,*)'do you want change lambda? y/[n]'
c	call flush(6)
c	read(*,'(a)')yesno
c	yesno='n'
c
c	go to read new lambda
c
	if(yesno.eq.'y')go to 1001

	write(*,*)'do you want plot log D.E.M.? y/n'
	call flush(6)
	read(*,'(a)')yesno
	call flush(5)
c
	write(4,*)' new value in mesh points:'
	write(4,*)(xt(i),i=1,nnodi)
	write(4,*)(hcf(i),i=1,nnodi)
c
c	si chiamma ttmax per sapere il max del prodotto g(T)*fi(T)
c
	do line=1,nn
	tttmax(line)=ttmax()
	demm(line)=demmax
	end do

c
c	si preparano con spline le DEM da plottare
c

	atfin=xt(nnodi)
	atin=xt(1)
	dpassi=(atfin-atin)/100.
	ddpassi=(atfin-atin)/1000.
	k=1
	xint(1)=atin
	xxint(1)=atin
333	do k=2,100
	xint(k)=xint(k-1)+dpassi
	end do
	do k=2,1000
	xxint(k)=xxint(k-1)+ddpassi
	end do
	do i=1,1000
	xx=xxint(i)
	call spline(xt,hcf,nnodi,2e30,2.e30,y2)
	call splint(xt,hcf,y2,nnodi,xx,y)
	write(99,*)xx,y
	end do
	do i=1,100
	xx=xint(i)
	call spline(xt,hcf,nnodi,2e30,2.e30,y2)
	call splint(xt,hcf,y2,nnodi,xx,y)
	yteor(i)= y
	end do
c	write(title,3333)rms
3333	format(' ~RMS=',f6.3)

c	write(*,*)'do you want plot log D.E.M.? y/n'
c	call flush(6)
c	read(*,'(a)')yesno

	if(yesno.eq.'y')then
	write(*,*) (xint(it),it=1,100)
	write(*,*) (yteor(it),it=1,100)
	call flush(6)
	end if

	write(*,*)'would you change lambda?'
	call flush(6)
	read(*,*)yesno
	if (yesno.eq.'y')then
	nnpass=1
	go to 1001
	end if


	RETURN

	write(*,*)'do you want plot OSS vr. Exp:y/[n]'
	read(*,110)yesno
	if(yesno.eq.'y')then
c	call plottaoss(dev1,yoss,yint,nn,'Obs'
c     +			,' Exp  ',title)
	read(*,'(a)')yesno
	write(*,*)'do you want plot on laser? y/n'
	read(*,110)yesno
	if(yesno.eq.'y')then
c	call plottaoss('ps',yoss,yint,nn,'Obs'
c     +			,' Exp  ',title)
	end if
	end if
	write(*,*)'do you want plot quality factor vr. Tmax: y/[n]'
	read(*,110)yesno
	if(yesno.eq.'y')then
	write(*,*)'enter xmin,xmax,ymin,ymax'
	read(*,*)x1,x2,y1,yy2
c	call plottaqf(dev1,tttmax,rip,nn,'log Tmax'
c     +			,' Ie/Io  ',title,x1,x2,y1,yy2,verr,yoss)
	read(*,'(a)')yesno
	write(*,*)'do you want plot on laser? y/n'
	read(*,110)yesno
	if(yesno.eq.'y')then
c	call plottaqf('ps',tttmax,rip,nn,'log Tmax'
c     +			,' Ie/Io  ',title,x1,x2,y1,yy2,verr,yoss)
	end if
	end if
	write(*,*)'do you want plot quality factor vr.catalog n:y/[n]'
	read(*,110)yesno
	if(yesno.eq.'y')then
	write(*,*)'enter xmin,xmax,ymin,ymax'
	read(*,*)x1,x2,y1,yy2
c	call plottaqf(dev1,sert,rip,nn,'catalog number'
c     +			,' Ie/Io  ',title,x1,x2,y1,yy2,verr,yoss)
cc	call plottaqf(dev1,all,rip,nn,'wavelenght'
c     +			,' Ie/Io  ',title,x1,x2,y1,yy2,verr,yoss)
	read(*,'(a)')yesno
	write(*,*)'do you want plot on laser? y/n'
	read(*,110)yesno
	if(yesno.eq.'y')then
c	call plottaqf('ps',sert,rip,nn,'catalog number'
c     +			,' Ie/Io  ',title,x1,x2,y1,yy2,verr,yoss)
	end if
	end if
c	write(66,*)nnodi
	do j=1,nnodi
c	write(66,*)xt(j),hcf(j)
	end do
	return
	end
	
	function fct()
c
c	calcolo del segnale teorico
	dimension y2(100),fi(100),pr(100),pr1(100)
	common/tt/ain0,ddt,mmt,att(100)
	common/spettro/al(100),asp(100,200),at(200)
	common/pinco/xxa(100,100),yya(100,100),nt(100)
	common/caio/yfi(100,100)
	common/line1/line,kkt
	common/cfff/nnodi,hcf(100),xt(100)
	common/ngauss/ngauss
	external func
	ain=ain0
	segnale=0
	nnn=nt(line)
	ain=xxa(2,line)
	afin=xxa(nnn-1,line)
c	write(41,*)'xt= ',(xt(i),i=1,nnodi)
	do jt=1,nnn
c	do jt=2,nnn-1
	yfi(jt,line)=0.
	xb=xxa(jt,line)
c	if(xb.gt.xt(nnodi))go to 10
	call spline(xt,hcf,nnodi,2.e30,2.e30,y2)
	call splint(xt,hcf,y2,nnodi,xb,fy)
c
c	log [G(T)*f(T)]
c
	fi(jt)=fy
	yfi(jt,line)=yya(jt,line)+fy
	pr(jt)=yya(jt,line)
	pr1(jt)=yfi(jt,line)
c	write(41,*)xb,jt,fy,pr(jt),pr1(jt)
c	write(*,*)xb,jt,fy,pr(jt),pr1(jt)
	end do
10	aain=10.**ain
	aafin=10.**afin
c	kpassi=20
c	call integra(func,aain,aafin,segnale,kpassi)
c	go to 541
	deltat=(afin-ain)/ngauss
	t1=ain
	do n=1,ngauss
	t2=t1+deltat
	aain=10.**t1
	aafin=10.**t2
	call qgaus(func,aain,aafin,ss)
	segnale=segnale+ss
	t1=t2
	end do
541	fct=segnale
	return
	end

c	serve ad ordinare in ordine crescente i valori dei mesh points
c
	subroutine ordina(xt,hfc,nnodi)
	dimension xt(100),hfc(100),xo(100),hfco(100),index(100)
	do i=1,nnodi
	xo(i)=xt(i)
	hfco(i)=hfc(i)
	end do
	call kb07a(xo,nnodi,index)
	do i=1,nnodi
	xt(i)=xo(i)
	hfc(i)=hfco(index(i))
	end do
	return
	end



