C	Inversione di n sistema con il metodo della Max entropia
c	aggiunto termine di smorzamento
c	si usa la spline
c
c	gennaio 1991
c 	aggiornato per XVision gennaio 1992
c
c	febbraio 1994 adattato per IDL
c
	dimension tlog(50),add(100)
	dimension de0(100),abb0(50),abb(50),el0(50)
	dimension all(100),verr(100),err(100),segnale(100)
      dimension  rde(100),kkk(100),ymed(100),ysum(100)
      dimension  de(100), dd(100),t1(100),t2(100),tttmax(100)
	dimension t_g_dem(100)
	dimension aatt(200),ssp(100,200),tmax(100),ll(100)
	dimension nsert(100),nnsert(100)
	external func
	external handler
	common/rip/rip(100)
	common/tt1/tt1
	common/nsert/nsert(100)
	common/line1/line,kkt
	common/cfff/nnodi,hcf(100),xt(100)
	common/spettro/al(100),asp(100,200),at(200)
	common/pinco/xxa(100,100),yya(100,100),nt(100)
	common/tt/ain0,ddt,mmt,att(100)
	common/ion/ela(100),ion(100),tmax,tttmax,demm(100)
	common/passo/passo
	common chi, sd1,sd
	common/ilog/ilog
	common/dev1/dev1
	common/ngauss/ngauss
	common/rd/rstar,dstar
	character*80 title
	character el1*2
	character el*5
	character el0*5
	character ela*5
	character aion*5
	character ion*5
	character dev1*5
	character dev2*5
	character flusso*80
	character abbund*80
	character spettro*80
	character inputdem*80
c	character emm*80
	character*1 yesno
	ierr= ieee_handler('set','all',handler)
	ierr= ieee_handler('clear','inexact',handler)
	nin=5
      	nout=6
	ngauss=4
	isert=1
	write(*,*) 'entra nome file flusso'
	call flush(6)
	read(*, '(a)')flusso
c	write(*,*) 'entra e i dati sono del SERTS(arcsec-2):yes=0'
c	read(*,*)isert
	isert=1
	write(*,*) 'entra nome file spettro'

	call flush(6)
	read(*, '(a)') spettro
c	abbund='aaa1.dat'
c	write(*,*) 'entra nome file abbondanze'
c	call flush(6)
c	read(*, '(a)')abbund
2225	format(a80)
	abbund='AAA.DAT'
      open(unit=2, file=flusso, status='old') 
	open(unit=33,file='AAA.DAT',status='old')
      open(unit=22, file=abbund, status='old') 
      open(unit=44, file=spettro, status='old') 
      open(unit=4, file='ris.dat') 
	do i=1,50
	abb0(i)=0.
	end do
	do i=1,30
	abb(i)=0.
	end do
	write(4,*)' output programma ENTROPIA '
	write(4,2225)'file obs = ',flusso,'file spp = ',spettro
	l=0
1112      l=l+1
c	read(2,*,end=111,err=111)el,aion,all(l),de(l),err(l),de0(l),tmax(l)
	read(2,*,end=111,err=111)nsert(l),el,aion,all(l),de(l),err(l),
     +	de0(l),tmax(l)
c	write(31,*)nsert(l),el,aion,all(l),de(l),err(l),
c     +  de0(l),tmax(l)

	if(isert.eq.0)de(l)=de(l)/1.9e-12
	if(isert.eq.0)verr(l)=err(l)/1.9e-12
	write(9,*)nsert(l),el,aion,all(l),de(l),err(l),de0(l),tmax(l)
	verr(l)=err(l)
	go to 1112
111	continue
	lll=l-1
	write(9,*)'lll = ',lll
c
c	legge abbondanze

c
	ja=0
2114	 read(33,1,end=2113,err=2113)el1,j,abb0(j)
1	format(a2,i3,f5.2)
	el0(j)=el1
c	write(*,*)j,el0(j),abb0(j)
	ja=ja+1
	go to 2114
2113	continue
	close(33)
	ja=0
2112	 read(22,*,end=2111,err=2111)j,el0(j),abb(j)
c	write(*,*)j,el0(j),abb(j)
	ja=ja+1
	go to 2112
2111	continue
	close(22)
	do j=1,30
	ab1=abb0(j)+abb(j)
	abfe=ab1-abb0(26)
	write(9,*)j,el0(j),abb0(j),abb(j),ab1,abfe
	write(4,*)j,el0(j),abb0(j),abb(j),ab1,abfe
	end do
	isegnale=0
c	write(*,*)'vuoi il segnale generato con errore??si=0'
c	read(*,*)isegnale
	if(isegnale.ne.0)then
	do 444 i=1,l
	de(i)=de0(i)
444	continue 
	write(4,*)'segnale costruito senza errore'
	end if
	read(44,*)t,dt,kkt,ei
555	format(f10.2,5(10f6.2))
	tt1=t
	i=0
	yesno='y'
c	write(*,*)'i flussi sono per sterad? [y]/n'
c	read(*,'(a)')yesno
	yesno='n'
2222	i=i+1
	read(44,*,end=222,err=222)nnsert(i),el,aion,al(i),(ssp(i,it),it=1,kkt)
	do jj=1,30
	if(el.eq.el0(jj))then
	continue
	do it=1,kkt
	ssp(i,it)=ssp(i,it)-abb(jj)
	end do
	endif
	end do
c
c	si divide per 4*pi per tener conto che i dati sono per sterad
c
	yesno='n'
	if(yesno.ne.'n')then
c
c!!!!	per i dati del sert si deve aggiungere 10.61976	, cioe' si
c	tiene conto che il flusso e' per arcsec**2
c
	do it=1,kkt
	if(isert.eq.1)ssp(i,it)=ssp(i,it)+1.099210
	if(isert.eq.0)ssp(i,it)=ssp(i,it)+1.099210+10.6179
	end do
	end if
	ion(i)=aion
	ela(i)=el
c	write(*,*)el,aion,al(i),(ssp(i,it),it=1,kkt)
	go to 2222
222	continue
	nn=i-1
C	
C	si voglion escludere alcune righe
C
	i=0
	yesno=' '
	write(*,*)' do you select lines from file?  [n]/y'
	call flush(6)
	read(*,'(a)') yesno
	if(yesno.eq.'y'.or.yesno.eq.'Y') then
c	write(*,*)'le lambda delle righe da escludere sono su file:y/[n]'
c	read(*,*)yesno
	yesno='y'
6532	if(yesno.eq.'y')then
	read(55,*,end=6531,err=6531)mmline
	write(9,*)mmline
	if(mmline.lt.0.)go to 44444
	do line=1,nn
	if(mmline.eq.nnsert(line))then
	i=i+1
	ll(i)=line
	end if
	end do
	go to 6532
	endif
55555	if(yesno.ne.'y')then
	write(*,*)' entra lambda riga da escludere'
	read(*,*)mmline
	if(mmline.lt.0.)go to 44444
	do line=1,nn
	if(mmline.eq.nnsert(line))then
	i=i+1
	ll(i)=line
	end if
	end do
	go to 55555
	endif
c		do jj=1,i
6531	continue
44444	kk=i
	write(9,*)kk,nn,(ll(k),k=1,nn)
	kline=0
		 do line=1,nn
	kline=kline+1
	ela(kline)=ela(line)
	de(kline)=de(line)
	dd(kline)=dd(line)
	all(kline)=all(line)
	t1(kline)=t1(line)
	t2(kline)=t2(line)
	tmax(kline)=tmax(line)
	err(kline)=err(line)
	verr(kline)=verr(line)
	al(kline)=al(line)
	ion(kline)=ion(line)
	nsert(kline)=nsert(line)

	do it=1,kkt
	ssp(kline,it)=ssp(line,it)
	end do

	do i=1,kk
	if(ll(i).eq.line)then
	kline=kline-1
	go to 12345
	end if

	end do
12345	continue
		end do
c	if(yesno.eq.'y')go to 6532
c	if(yesno.ne.'y') go to 55555
c	go to 55555
	nn=nn-kk
	write(*,*)nn,kkt
	call flush(6)
	end if
	write(9,*)kline,nn,(ll(k),k=1,nn)
	write(title,'(a)')el
c
c PLOT G(T)
c
c	write(*,*) ' quale device grafico vuoi? (x: o tt:)'
c	read(*,'(a)') dev1
	dev1='x:'
	dev2='ps:  '
	write(*,*) ' do you wish to plot G(T)? [N]/Y'
	call flush(6)
	read(*,'(a)') yesno


c	if(yesno.eq.'y'.or.yesno.eq.'Y') then

c	if(yesno.eq.'y'.or.yesno.eq.'Y') then
cc
cc
cc	T vr. ssp  creare il vettore T
cc
	tlog(1)=tt1+6
	do it=2,kkt
	tlog(it)=tlog(it-1)+dt
	end do
	
	write(6,*) (tlog(it),it=1,kkt)
	call flush (6)
	do it=1,kkt
	write(6,*) (ssp(i,it),i=1,nn)
	call flush(6)
	end do
cc
cc
c	endif
c	endif
c
c
c PLOT SEGN/G(T)
c
	do line=1,nn
	rde(line)=de(line)
	enddo
c
	write(*,*)'enter Rstar(rsun), Dstar(pc)'
	callflush(6)
	read(*,*)rstar,dstar

	rstar=rstar*6.9e10
	dstar=dstar*3.08e18
	write(*,*) ' do you wish to plot SIGN/G(T)? [n]/y'
	call flush(6)
	read(*,'(a)') yesno
c	if(yesno.eq.'y'.or.yesno.eq.'Y') then

	write(6,*) (rde(i),i=1,nn)
	call flush(6)

c	call plottaEM(dev1,ssp,nn,kkt,rde,dt)
c	read(*,'(a)') yesno
c	write(*,*) 'vuoi copia sulla laser? [N]/Y'
c	read(*,'(a)') yesno
c	if(yesno.eq.'y'.or.yesno.eq.'Y') then
c	call plottaEM('ps',ssp,nn,kkt,rde,dt)
c	endif
c	endif


c
c
c	nn= numero dati ingresso
c
c
	iit=kkt
	do 4440 i=1,nn
	do 4441 kt=1,iit
	asp(i,kt)=ssp(i,kt)
4441	continue 
4440	continue	
	at(1)=t
	if(t.lt.0)at(1)=tt1+6
	do  9 it=1,iit-1
	at(it+1)=at(it)+dt
9	continue
c
c
98765	continue
c	write(*,*)'log TEMP INIZ.'
c	read(*,*)ain0
	ain0=4.3
	ddt=1.
	passo=1.
c
	write(*,*)'enter iit= # steps of T,  % ERROR if < 0 gets 
     &	the observed value'
	call flush(6)
	read(*,*)iit,error
	if(error.ge.0)then
	do 4442 line=1,nn
	verr(line)=error*de(line)
c	err(line)=error*de(line)
4442	continue	
	end if
	err1=0
	do line=1,nn
	err1=err1+verr(line)
c	err1=err1+err(line)
	end do
	errm=err1/nn
	go to 99991
	do line=1,nn
	if(verr(line).lt.0.01*errm)then
c	if(err(line).lt.0.01*errm)then
	verr(line)=errm*.01
c	err(line)=errm*.01
	end if
	end do
99991	continue
	afin0=at(iit)
	mmt=.01+(afin0-ain0)/ddt
	write(9,*)at(iit),at(1),ddt,mmt,iit
	write(4,*)at(iit),at(1),iit

c
c	si calcola integrale g(T) [ fi(T)=1sempre]
c
	nnodi=9
	xt(1)=at(1)
	do n=2,nnodi
	xt(n)=xt(n-1)+0.5
	end do
	do i=1,nnodi
	hcf(i)=0.
	end do


	do it=1,iit
	ymed(it)=0.
	ysum(it)=0.
	kkk(it)=0
	end do

c	
c	entra se si vuole che intervallo integrazione sia scelto
c	in modo che sia n ordini sotto il picco della gt
c
c	write(*,*)'integration limit: 0[logT=4 to 8]'
      	write(*,*)' n [orders of magnetudine less then  
     &	peak value of GT]'
	call flush(6)
	read(*,*)itest
	do 91 line=1,nn
	newgt=funx(itest)
	ain=ain0
	x0=line
	dd(line)=fct()
	nnn=nt(line)
	ai=xxa(2,line)
	af=xxa(nnn-1,line)
	t1(line)=xxa(2,line)
	t2(line)=xxa(nnn-1,line)
	dddd=alog10(de(line)/dd(line))
	write(9,*)all(line),de(line),dd(line),dddd,ai,af,
     +	tmax(line),nsert(line)
c	write(*,*)de(line)/dd(line),tmax(line)
c	write(31,*)ela(line),ion(line),all(line),de(line),dd(line),
c     +	de(line)/dd(line),tmax(line),t1(line),t2(line)
	do it=1,iit
	if(tmax(line).ge.at(it).and.tmax(line).le.at(it+1))then
	ysum(it)=ysum(it)+de(line)/dd(line)
	kkk(it)=kkk(it)+1
	end if
	end do

91 	continue

	nll=0
	do it=1,iit
	if(ysum(it).gt.0)then
	ymed(it)=log10(ysum(it)/kkk(it))
	write(9,*)at(it),ymed(it)
	nll=nll+1
	endif
	end do
	write(9,*)nll
c	write(17,*)nll
	do it=1,iit
	if(ysum(it).gt.0)then
c	write(17,*)at(it),ymed(it)
	endif
	end do
	close (17)

1111	format(10(1pe10.3))

c	
c	PLOT SEGN/INTEGRALE vr. Tmax
c
	write(*,*)'do you wish to plot SIGN/G(T) vr. Tmax? [n]/y'
	call flush(6)
	read(*,'(a)') yesno

	write(6,*)(tmax(i),i=1,nn)
	write(6,*)(dd(i),i=1,nn)
	call flush(6)



22222	continue
c	write(6,*)'enter mesh point from file? [y]/n'
c	read(*,'(a)') yesno
	yesno='y'
	do kk=1,100
	xt(kk)=0.
	end do
	if(yesno.eq.'y')then
	write(6,*) 'entra nome file inputdem'
	call flush(6)
	read(*, '(a)') inputdem
      open(unit=15, file=inputdem, status='old') 
	write(4,*)inputdem
	write(9,*)inputdem
	read(15,*)nnodi
	do i=1,nnodi
	read(15,*)xt(i),hcf(i)
	end do
	rewind 15
	close (15)
	go to 6666
	end if
c
c	entrano numero mesh points e coordinate
c
	write(*,*)'enter #mesh points '
	read(*,*)nnodi
	write(*,*)'enter values T,f(T)'
	write(16,*)nnodi
	do i=1,nnodi
	read(*,*)xt(i),hcf(i)
	write(16,*)xt(i),hcf(i)
	end do
	rewind 16
	close (16)


6666	continue

c
c
c	si chiama la subroutine per l'inversione
c
	write(9,*)(hcf(i),i=1,nnodi)
	call inv_entr(all,de,nn,verr,verr)
	do line=1,nn
        dd(line)=fct()
        t_g_dem(line)=fff()
        teff=alog10(t_g_dem(line)/dd(line))
        dem_t=demeff(teff)



	add(line)=alog10(de(line)/dd(line))
	write(31,4321),line,ela(line),ion(line),all(line),de(line),dd(line),
     +	add(line),tmax(line),t1(line),t2(line),
     +	tttmax(line),demm(line),rip(line),verr(line), teff,dem_t

4321	format(i4,2(a5),1x,f6.2,2(e11.3),7(f6.2),e11.3)
	end do
	STOP

c	call inv_entr(all,de,nn,err,verr)
c	
c	PLOT SEGN/INTEGRALE vr. Tmax del prodotto g(T)*fi(T)
c
c	write(*,*) ' vuoi plot SIGN/G(T) vr. TTTmax? [N]/Y'
c	read(*,'(a)') yesno
c	if(yesno.eq.'y'.or.yesno.eq.'Y') then
c	call plottaTM(dev1,tttmax,t1,t2,dd,nn,'log T','E.M.',title,de)
c	read(*,'(a)') yesno
c	write(*,*) 'vuoi copia sulla laser? [N]/Y'
c	read(*,'(a)') yesno
c	if(yesno.eq.'y'.or.yesno.eq.'Y') then
c	call plottaTM('ps',tttmax,t1,t2,dd,nn,'log T','E.M.',title,de)
c	endif
c	endif
	do line=1,nn
	add(line)=alog10(de(line)/dd(line))
	write(31,4321)ela(line),ion(line),all(line),de(line),dd(line),
     +	add,tmax(line),t1(line),t2(line),
     +	tttmax(line),demm(line),rip(line),verr(line)
c4321	format(2(a5),1x,f6.2,2(e11.3),7(f6.2),e11.3)
	end do
	write(*,*)(tttmax(iline),iline=1,nn)
	write(*,*)(add(iline),iline=1,nn)
	call flush(6)
	STOP
	write(*,*)'do you wish to change mesh points?'
	read(*,'(a)') yesno
	if(yesno.eq.'y')go to 22222

	stop
	end

	integer function handler(sig,code,contex)
	integer sig,code,contex(5)
	write(*,*)'ci sono!!'
c	call abort()
	return
	end 


