;+ ***********************************************************************
; NAME:
;	PROFDECOMP1D
;
; PURPOSE:
;  PROCEDURE DE TRACE D'UN PROFIL RH UNIDIMENSIONNEL ET DE SA DECOMPOSITION 
;  EN SOURCES ELEMENTAIRES (INTENSITE OU POLARISATION) A UN TEMPS DONNE
;  A PARTIR D'UN TABLEAU CONTENANT DES IMAGES EN FONCTION DU TEMPS PROVENANT
;  DE LA FONCTION LEC1D (IMAGE 1D DU RH) (TABLEAUX DE 2D "INTENSITE/TEMPS" OU
;  3D "INTENSITE+POLARISATION/TEMPS"). 
;
; CATEGORY:
;	NRH1 Positions
;
; CALLING SEQUENCE:
;	PROFDECOMP1D, Tabl_X, Lun, Gliss, Poids, Npos
;
; INPUTS:
;      TABL_X ..... Tableau issu de LEC1D contenant les images(2 dimensions) 
;                   tabl_X(num_Canal,Temps)(3 Dim sinon )
;                           		tabl_X(Num_Canal,Temps,I:V)
;      LUN .......  unite logique assignee au fichier
;	 GLISS ...... Nombre de canaux sur lequel la moyenne
;					glissante est calculee
;
;	 POIDS ...... Poids minimum d'une source a retenir (->
;					SIGRMN)
;        NPOS ....... Nombre maximum d'iterations sur 1 image
;	
; KEYWORD PARAMETERS:
;       HEURE ...... heure du profil à tracer,
;                                   format : chaine de caractere hh:mm:ss:ms
; 
;                                    si on est en mode animation
;       CANAUX ..... Tableau une dimension, deux elelments
;				   contient les canaux de 
;				   debut et de fin de visualisation
;	NOERASE .... Permet de forcer le passage au mode 
;				   NOERASE si on veut visualiser un ou
;				   plusieurs profils ne se suivants pas 
;				   dans la meme fenetre
;	WIN ........ Specifie le numero de fenetre dans lequel a
;
;       LASER ...... Permet de rediriger la sortie graphique sur
;                                   l'imprimante laser
;       INTENSITE .. Trace d'un profil d''intensite
;       POLAR........Trace d'un profil polarise
;       SIMPLE ..... Ajustement d'une Gaussienne
;       CONST, DROITE, PARABOLE ... Gaussienn + cst., droite 
;
; OUTPUTS:
;	Trace
;
; COMMON BLOCKS:
;	NRH!_GLOBAL
;
; PROCEDURE:
;  METHODE: CALCUL D'UNE MOYENNE GLISSANTE SUR gliss CANAUX, IDENTIFICATION
;  DES POINTS QUI DEPASSENT LA MOYENNE GLISSANTE (diff image - image lissee,
;  normalisee sur 1, >= poids). AJUSTEMENT D'UNE GAUSSIENNE (OU D'UNE SOMME
;  D'UNE GAUSSIENNE ET D'UN POLYNOME D'ORDRE 0, 1 OU 2) A LA SOURCE AINSI 
;  IDENTIFIEE.
;
; EXAMPLE:
;	Please provide a simple example here
;
; MODIFICATION HISTORY: (bonmartin@obspm.fr)
;	16/11/98	adapte du logiciel XHELIO
;	09/04/99	ajout mot cle message: liste des pos. identifiees 
;       31/08/00	ajout possibilite d'ecrire resultats de la decompo-
;			sition en sources elementaires dans fichier ASCII
;			(nom standard: 'nrhascipos1d.dat')
;-*******************************************************************



PRO profdecomp1d,tabl_X, lun, gliss, poids, npos, laser=Laser, Tek=Tek $
             ,Intensite=intensite, Polar=Polar, WIN=win, NOERASE=noerase $
             ,BORD=bord,HEURE=heure, CANAUX=canaux $
             ,FILTER=filter, BANDE=bande, SIMPLE=simple, CONST=const,$
		DROITE=droite, PARABOLE = parabole, MESSAGE=message, $
		WRITEASCI=writeasci

;* ---------------------------- *
;*	Test de la syntaxe      *   à revoir !!!!!!
;* ---------------------------- *
;
IF N_PARAMS() LE 1 THEN BEGIN
   print,$
    'Syntaxe: profdecomp1d, Tableau a tracer, Numero d unite logique'$
   , '          [,LASER=laser, TEK=tek,' $ 
   , '            ,Intensite=intensite,Polar=Polar, WIN=win,NOERASE=noerase,' $
   , '            ,BORD=bord,HEURE=heure, CANAUX=canaux,' $
   , '		]'
; 
ENDIF

;
;* ------------------------------------------- *
;*	Declaration des variables globales     *
;* ------------------------------------------- *
;            

COMMON NRH1_GLOBAL,nomfich,FichInfo      

;
; Memorisation des informations principales
; sur le fichier
;
file = fichinfo(lun)
heuredeb = file.Heurdeb
heurefin = file.HeurFin
dt = file.paspix(1)	; intervalle de temps
dx = file.paspix(0)	; intervalle de canaux
cdeb = file.valorg(0)
cfin = cdeb + (file.nbpix(0)-1)*dx

;
; --------------------------------------------
;       Determination des valeurs par defaut
; --------------------------------------------
IF (N_ELEMENTS(gliss) EQ 0) THEN gliss = 7
IF (N_ELEMENTS(poids) EQ 0) THEN poids = 0.5 
;
; ---------------------------------- 
;	Transformation de l'heure     
;         en milliseconde            
;- ---------------------------------- 
;

IF NOT KEYWORD_SET(heure) THEN h = heuredeb ELSE h = HMS(heure)

;---------------------------------------------
; Test sur la conformite de l'heure specifiee
;---------------------------------------------

IF (h LT heuredeb OR h GT heurefin) THEN BEGIN
	Print,[ 'Heure non conforme : elle doit etre comprise entre ', $
           MSH(heuredeb),' et ',MSH(heurefin), '.',$
           'Je prends Heure de depart = ',MSH(Heuredeb)]  
	h = heuredeb

ENDIF 

;-------------------------------------------
; Test sur la conformite des canaux specifies
;-------------------------------------------

IF NOT KEYWORD_SET(canaux) THEN BEGIN
	cnldeb=cdeb
	cnlfin=cfin
ENDIF ELSE BEGIN
	cnldeb=canaux(0)
	cnlfin=canaux(1)
ENDELSE
IF (cnldeb GE cnlfin) THEN BEGIN
	print,'CANAL DEB > OR = THAN CANAL FIN'
	RETURN
ENDIF
IF(cnldeb GT cdeb AND cnldeb LT cfin) THEN cnldeb=cdeb+FIX((cnldeb-cdeb)/dx)*dx
IF(cnlfin GT cdeb AND cnlfin LT cfin) THEN cnlfin=cfin+FIX((cnlfin-cfin)/dx)*dx
canaux = [cnldeb, cnlfin]

xcoor = cdeb + findgen(file.nbpix(0))*dx	; Tableau d'Abscisses: canaux

num = indice(heuredeb,h,dt, heurelue)
heure = MSH(heurelue)
h = heurelue

;---------------------------------------------------------------
;  CHARGEMENT DES INFORMATIONS A AFFICHER
;---------------------------------------------------------------

freq=motcle(lun,'FREQ')
label=motcle(lun,'LABEL')
type=motcle(lun,'TYPE')
xtitle = 'Canal /' + string(label)
ytitle = 'Amplitude /' + string(freq)
subtitle = string(fichinfo(lun).date) + '   ' + heure

;-------------------------------------
; Test de la specification d'un numero
; de fenetre
;-------------------------------------
IF (NOT KEYWORD_SET(win)) THEN BEGIN
  win = intarr(2)
  win(0) = 1
  win(1) = 2
ENDIF


  ;formatage de la fenetre de visualisation

   ;* ----------------------------------------------- *
   ;*	Prise en compte du périphérique de sortie     *
   ;* ----------------------------------------------- *

IF KEYWORD_SET(laser) THEN nrh1_postscript,/ON	


     ;------------------------------------
     ;    Selection du profil
     ;------------------------------------  
     
IF (!P.MULTI(1) NE 0) THEN Taille=0.4 ELSE Taille=0.8


;	IF (KEYWORD_SET(bord)) THEN $
;	IF NOT KEYWORD_SET(LASER) THEN $
;	WINDOW ,win(0),title='INTENSITE      Fenetre numero :'+string(win(0)),$
;            retain=2,xpos=520,ypos=350,xsize=500,ysize=500

IF(KEYWORD_SET(intensite)) THEN BEGIN
	LoadIntens,Tabl_x,imagex,num,Fichinfo,lun
	ytitle = 'I - ' + ytitle
ENDIF ELSE IF(KEYWORD_SET(polar)) THEN BEGIN
	LoadPolar, Tabl_x,imageX,num,FichInfo,lun,ok
	IF (ok NE 1) THEN BEGIN
		print, "*** Le fichier ne contient pas d'image V ***"
		return
	ENDIF
	ytitle = 'V - ' + ytitle
ENDIF 
imagex1 = imagex

maximage = MAX(imagex, MIN = minimage)
space = 0.1*(maximage-minimage)
yrange = [minimage-space, maximage+space]

red = [0, 1, 1, 0, 0, 1] & green = [0,1,0,1,0,1] & blue = [0,1, 0, 0, 1, 0]
TVLCT, 255*red, 255*green, 255*blue

;----------------------------------------------
; Trace du profil brut (en blanc sur l'ecran) :
;----------------------------------------------
ymargin = [20,2]	; revoir ce choix ***
PLOT,xcoor,imagex, $
		XTITLE = xtitle, YTITLE = ytitle, SUBTITLE = subtitle, $
		xstyle=1 ,ystyle=1, XRANGE=canaux, YRANGE=yrange,$
		YMARGIN=ymargin, NOERASE=noerase, COLOR=1, CHARSIZE=Taille
;----------------------------------------------
; Decomposition en sources Gaussiennes
;----------------------------------------------

decomp1d, xcoor, imagex, gliss, poids, POLAR=polar, $
	SIMPLE=simple, CONST=const, DROITE=droite, PARABOLE = parabole, $
	fitparam, dfitparam, sumgauss, npos


IF KEYWORD_SET(SIMPLE) THEN text1 = ' (Gaussian)'
IF KEYWORD_SET(CONST) THEN text1 = ' (Gaussian + constant)'
IF KEYWORD_SET(DROITE) THEN text1 = ' (Gaussian + straight line)'
IF KEYWORD_SET(PARABOLE) THEN text1 = ' (Gaussian + parabola)'
message=STRARR(npos>1)
IF npos EQ 0 THEN  message(0) = 'No local source found '

chantxt = STRARR(10) & fwhmtxt = STRARR(10) & amplitxt = STRARR(10)
fluxtxt = STRARR(10) & intervtxt = STRARR(10)

FOR ipos = 0, npos-1 DO BEGIN
	ampli = fitparam(0, ipos)
	chan = fitparam(1, ipos)
	fwhm = fitparam(2, ipos)
	flux = 1.06446*ampli*fwhm

;	writeasci=1
	IF KEYWORD_SET(writeasci) THEN BEGIN
		fichier='nrhascipos1d.dat'
		a = FINDFILE(fichier, COUNT=new)
		IF new EQ 0 THEN BEGIN
        		openw,1,fichier, ERROR=err 
			PRINTF, 1, FORMAT='(4X, "Array", 1X, "Frequency", 11X, "UT", 1X, "Chann.", 3X, "FWHM", 4X, "Amplitude", 1X, "Flux density", /)'
		ENDIF ELSE BEGIN
        		openu, 1, fichier, ERROR=err
			tmp = FSTAT(1)
			POINT_LUN, 1, tmp.size ; positionner pointeur a la fin du fichier
		ENDELSE
		IF err NE 0 THEN RETURN
		text= STRING(label) + ' ' + STRING(freq)  + ' ' + heure + $
			STRING(FORMAT = '(F7.2)', chan)  + $
			STRING(FORMAT = '(F7.2)',fwhm) + $ + 
			STRING(ampli) + STRING(flux)
		PRINTF, 1, text
		CLOSE, 1
	ENDIF


; Impression des resultats sur l'ecran
	message(ipos) = STRCOMPRESS("Channel: " + STRING(chan)+ $
		"  Width: " + STRING(fwhm) + $
		"  Amp: " + STRING( ampli) + " Flux: " + STRING(flux) )

	errflux = flux*sqrt( (dfitparam(0, ipos)/ampli)^2 $
		+ (dfitparam(2, ipos)/fwhm)^2 )
	chantxt(ipos) = STRING(FORMAT = $
		 '(F7.2, "/", F5.3, 1X)', chan, dfitparam(1, ipos))+'!C'
	fwhmtxt(ipos) = STRING(FORMAT= $
		'(F5.2,"/",F5.3)',fwhm, dfitparam(2, ipos))+'!C'
	amplitxt(ipos) = STRING(FORMAT= $
		'(1X, G9.3, "/", G9.3)', ampli, dfitparam(0, ipos))+'!C'
	fluxtxt(ipos) = STRING(FORMAT= $
		'(1X, G9.3, "/", G9.3)', flux, errflux)+'!C'

ENDFOR

;----------------------------------------------
; Inscription des resultats au-dessous du trace
;----------------------------------------------
text1 = STRING(FORMAT = '( "Sources with weight", F4.1)', poids)$
				 + STRING(FORMAT = $
	'(" above the running mean over ", F5.2)', gliss) $
	    + " channels" + text1
text1 = STRCOMPRESS(text1)
chars = FLOAT(!d.x_size) / FLOAT((STRLEN(text1)*!d.x_ch_size))
IF(chars GT 0.8) THEN chars = 0.8
annoty =  (ymargin(0)-10)*!d.y_ch_size
xyouts, 1,annoty, text1, /DEVICE, charsize = chars

IF npos GT 0 THEN BEGIN	;	Inscription si au moins une source localisee
	chantxt = STRCOMPRESS(chantxt) & chanlngth = MAX (STRLEN(chantxt))
	fwhmtxt = STRCOMPRESS(fwhmtxt) & fwhmlngth = MAX (STRLEN(fwhmtxt))
	amplitxt = STRCOMPRESS(amplitxt) & amplilngth = MAX (STRLEN(amplitxt))
	fluxtxt = STRCOMPRESS(fluxtxt) & fluxlngth = MAX (STRLEN(fluxtxt))
	intervtxt = STRCOMPRESS(intervtxt) & intervlngth = MAX(STRLEN(intervtxt))
	chars = FLOAT(!d.x_size)/(Float(chanlngth + fwhmlngth + amplilngth + $
		fluxlngth + intervlngth)*!d.x_ch_size)
	IF(chars GT 0.8) THEN chars = 0.8

	chtxt ='Channel!C' & fwtxt='FWHM!C' & amtxt = 'Amplitude!C'
	fltxt ='Flux density!C' & intxt = 'Interval!C'
	FOR I = 0, Ipos-1 DO BEGIN
		chtxt = chtxt + '!C' + chantxt(i)
		fwtxt = fwtxt +'!C'  + fwhmtxt(i)
		amtxt = amtxt +'!C'  + amplitxt(i)
		fltxt = fltxt +'!C'  + fluxtxt(i)
		intxt = intxt +'!C'  + intervtxt(i)
	ENDFOR
	annotx = 0 & annoty = annoty - 2*!d.y_ch_size * chars
	charsy = annoty / FLOAT(2*(npos+1)*!d.y_ch_size)
	IF chars GT charsy THEN chars = charsy
	xyouts, annotx, annoty, chtxt, /DEVICE, charsize=chars
	annotx = annotx + chars*chanlngth*!d.x_ch_size
	xyouts, annotx, annoty, fwtxt, /DEVICE, charsize=chars
	annotx = annotx + chars*fwhmlngth*!d.x_ch_size
	xyouts, annotx, annoty, amtxt, /DEVICE, charsize=chars
	annotx = annotx + chars*amplilngth*!d.x_ch_size
	xyouts, annotx, annoty, fltxt, /DEVICE, charsize=chars
	annotx = annotx + chars*fluxlngth*!d.x_ch_size
;	xyouts, annotx, annoty, intxt, /DEVICE, charsize=chars
ENDIF

;----------------------------------------------
; Traces de la somme des Gaussiennes (jaune) et
; de leur difference avec le profil brut (vert)
;----------------------------------------------
	if(ipos gt 0) then begin
		nuls = WHERE (ABS(sumgauss) LT 1.E-34)
		IF( (size(nuls))(0) GT 0) THEN sumgauss(nuls) = 0.
		IF KEYWORD_SET(laser) THEN BEGIN
			oplot, xcoor, sumgauss, THICK=3
			oplot, xcoor, imagex1 - sumgauss, LINESTYLE = 3
		ENDIF ELSE BEGIN
			oplot, xcoor, sumgauss, color = 5 
			oplot, xcoor, imagex1 - sumgauss, color = 3
		ENDELSE
	endif else begin
; Trace de la moyenne glissante si aucune source loacele n'a ete identifiee
		xyouts, .05, .1, '*** No local sources identified ***', $
			charsize = 0.6
		oplot, xcoor, mscan, color = 5
	endelse


IF KEYWORD_SET(laser) THEN BEGIN
    nrh1_postscript,/OFF
ENDIF
 
END