;      Procedure to plot output of programm DEMCDS.f
;
;       pro plotdem
;
spawn,  'demcds', unit=pipe, /noshell
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; enter name file fluxes
flusso=' '
readf, pipe, flusso
print,flusso
read,flusso
print,flusso
printf, pipe, flusso

;  enter name file G(T)
gtfile='  '
readf, pipe, gtfile
print, gtfile
read,gtfile
print,gtfile
printf, pipe, gtfile

;lines from file
enter="  '
yesno="  '
readf, pipe,  enter
print,enter
read,yesno
printf,pipe,yesno

;G(T) plot???
readf, pipe,nn,kkt
print, 'nn=',nn, 'kkt= ',kkt

readf, pipe,  enter
print,enter
read,yesno
printf,pipe,yesno

tlog=fltarr(kkt)
sp=fltarr(nn)
ssp=fltarr(nn,kkt)

print, 'nn= ',nn, 'kkt= ',kkt

readf, pipe, tlog

for it=0,kkt-1 do begin
readf, pipe, sp
	for i=0,nn-1 do begin
	ssp(i,it)=sp(i)
	endfor
endfor

;print, tlog, ssp

sp=fltarr(kkt)

;!p.background=255
!p.color=0
if (yesno eq 'y') then begin
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;PLOT G(T)

set_plot, 'x'
xmin=4
xmax=8
ymin=-17
ymax=-11

plot,tlog,sp, xtitle=' T',ytitle= ' F',$
xrange=[xmin,xmax],yrange=[ymin,ymax],$
 back=255, color=0 

	for i=0,nn-1 do begin
for it=0,kkt-1 do begin
	sp(it)=-ssp(i,it)
	endfor
ii=i
if(ii gt 5)then ii=ii-6
oplot,tlog,sp, linestyle=ii,color=0
endfor

y=' '
print, 'do you want to print on Laser? y/[n]'
read, yesno
if (yesno eq 'y') then pshard,file='gt.ps'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
endif

; enter Rstar and Dstar in pc
readf, pipe, enter
print,enter
read,rstar,dstar
print,rstar,dstar
printf, pipe, rstar,dstar

        rstar=rstar*6.9e10
        dstar=dstar*3.08e18
	ddrr=(dstar/rstar)*(dstar/rstar)

;Signal/G(T)????
readf, pipe,  enter
print,enter
read,yesno
printf,pipe,yesno

rde=fltarr(nn)
yy=fltarr(kkt)
readf, pipe, rde

if (yesno eq 'y') then begin
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;PLOT SIGNAL/G(T)

set_plot, 'x'
xmin=4
xmax=8
ymin=30
ymax=40

plot,tlog,yy, xtitle=' log T',ytitle= '  EM*T',$
xrange=[xmin,xmax],yrange=[ymin,ymax],$
back=255,color=0

        for i=0,nn-1 do begin
for it=0,kkt-1 do begin

yy(it)=alog10(rde(i))+ssp(i,it)+alog10(ddrr)+tlog(it)
        endfor
ii=i
if(ii gt 5)then ii=ii-6
oplot,tlog,yy, linestyle=ii,color=0
endfor

y=' '
print, 'do you want to print on Laser? y/[n]'
read, yesno
if (yesno eq 'y') then pshard,file='emt.ps'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
endif

;enter iit,error
readf, pipe,  enter
print,enter
read,iit,error
printf,pipe,iit,error

;enter # of orders of magitude
readf, pipe,  enter
print,enter
read,itest
printf,pipe,itest

;Signal/IntegralG(T)?????
readf, pipe,  enter
print,enter
read,yesno
printf,pipe,yesno

tmax=fltarr(nn)
dd=fltarr(nn)
readf, pipe,  tmax
readf,pipe,dd

if (yesno eq 'y') then begin
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;PLOT SIGNAL/INTEGRAL OF G(T)

set_plot, 'x'
xmin=4
xmax=8
ymin=19
ymax=24

plot,tmax,yy, xtitle=' log Tmax',ytitle= '  EM',$
xrange=[xmin,xmax],yrange=[ymin,ymax]

        for i=0,nn-1 do begin

yy(i)=alog10(rde(i)/dd(i))+alog10(ddrr)
oplot,tmax,yy, psym=1
endfor

y=' '
print, 'do you want to print on Laser? y/[n]'
read, yesno
if (yesno eq 'y') then pshard,file='emtmax.ps'

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
endif

;enter name of file of first DEM
inputdem=' '
readf, pipe,  enter
print,enter
read,inputdem
printf,pipe,inputdem

; si entra nella sub inv_entr
print, ' it starts inv_entr' 
;the mesh points are parameters?[y/n
readf, pipe,  enter
print,enter
read,yesno
printf,pipe,yesno

;normalize the value in the mesh points? y/[n]
readf, pipe,  enter
print,enter
read,yesno
printf,pipe,yesno

;;;;;;;;;;;;;;;;;;;;
;it computes the solution for lambda=0
;and stores it
;;;;;;;;;;;;;;;;;;

;enter lambda and #steps

;READ LAMDA AND # OF STEPS
readlamb:
readf, pipe,  enter
print,enter
read,alambd,nnpassi
printf,pipe,alambd,nnpassi


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;write chi2 after nn steps
readf, pipe, chi2,rms,alamd,bestchi
print, 'chi2= ', chi2,'R.M.S.= ',rms,' lambda= ',alamd,'best chi= ',bestchi

;ask for new lambda
readf, pipe,  enter
print,enter
read,yesno
if (yesno eq 'y') then begin
printf,pipe,yesno

goto, readlamb
endif

if (yesno ne 'y') then printf,pipe,yesno

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;output of the fit
nsert=intarr(nn)
ela=strarr(nn)
ion=strarr(nn)
all=fltarr(nn)
yoss=fltarr(nn)
yint=fltarr(nn)
verr=fltarr(nn)
per=fltarr(nn)
rip=fltarr(nn)
tmax=fltarr(nn)

readf,pipe,nsert
readf,pipe,format='((a))', ela
readf,pipe,format='((a))',ion
readf,pipe,all
readf,pipe,yoss
readf,pipe,yint
readf,pipe,verr
readf,pipe,per
readf,pipe,rip
readf,pipe,tmax

for l=0,nn-1 do begin
print,nsert(l),ela(l),ion(l),all(l),yoss(l),yint(l),verr(l),per(l),$
rip(l),tmax(l)
endfor


;' lambda = ',alamd,tttest,xttest
tttest=strarr(1)
xttest=strarr(1)

readf, pipe,format='(e10.4,2a)', alamb,tttest,xttest
print,'lambda= ',alamb,tttest,xttest



;ask for start again with the best solution
readf,pipe,enter
print,enter
read,yesno
printf,pipe,yesno
;
;do you want change lambda? y/[n
;readf,pipe,enter
;print,enter
;read,yesno
;print,yesno
yesno='n'

if (yesno eq 'y') then begin
printf,pipe,yesno
 goto, readlamb
endif


enter2='   '
yesno=' '
;if (yesno ne 'y') then begin
;  plot log D.E.M.
readf,pipe,enter2
print,enter2
read,yesno
printf,pipe,yesno
;
;enter dem(T)
xint=fltarr(100)
yteor=fltarr(100)
readf,pipe,xint
readf,pipe,yteor

yteor=yteor+alog10(ddrr)
set_plot, 'x'
xmin=4
xmax=8
ymin=20
ymax=25
plot,xint,yteor, xtitle='log T', ytitle=' log DEM',$
xrange=[xmin,xmax],yrange=[ymin,ymax]

oplot,xint,yteor, linestyle=0

;endif

;;;;;;;;;;;;;;;;;;

readf,pipe,enter
print,enter
read,yesno
if(yesno eq 'y') then begin
printf,pipe,yesno
goto, readlamb
endif

end