cd, 'k:\Granulation\doklad-pulkovo2008'
device, decomposed=0
datdir = DIALOG_PICKFILE(/DIRECTORY,TITLE='choose dir with dat files for wavelet analysis');,GET_PATH=new_dir)
CD, datdir
datfiles=FINDFILE('*i.dat',count=nfiles)
if nfiles EQ 0.0 then datfiles=FINDFILE('*V.DAT',count=nfiles)
 if nfiles eq '' then begin
       		print, 'no input files'
       		goto,kon
       		endif

datfiles=datfiles(sort(datfiles))
print, nfiles, '  files found: '
print, datfiles
;loadct,27


;---------------------------cycle starts-----------------------
for ii=0, nfiles-1 do begin
CLOSE,1
CLOSE,2
CLOSE,3

aa=read_ascii(datfiles(ii),count=N0)
if N0 eq 0.0 then begin
       		print, 'no input fields in data file:', datfiles(ii)
       		goto,kon1
       		endif

;xx=aa.field1(0,*)
;yy=aa.field1(1,*)

IF NOT FINITE(aa.field1(0,0)) THEN begin
PRINT, 'NaN in dat file found'
xx=aa.field1(0,1:N0-1)
yy=aa.field1(1,1:N0-1)
N0=N0-1
endif else begin
xx=aa.field1(0,*)
yy=aa.field1(1,*)
endelse

print, 'file: ', datfiles (ii), ' / number of points:  ',N0
print, xx(*),yy(*)
interv=xx-shift(xx,1)
shag=min(abs(xx-shift(xx,1)))
shagmax=max(abs(interv(1:n0-1)))

if shagmax ne shag then begin
	num=(xx(n0-1)-xx(0))/shag +1
	N0=num
	x1=findgen(n0)*shag+xx(0)
	y1=interpol(yy,xx,x1)
	;x=fltarr(n0+1)
	;y=fltarr(n0+1)
	x=dblarr(n0+1)
	y=dblarr(n0+1)
	x(1:n0)=x1(*)
	y(1:n0)=y1(*)
	window,0,xsize=600,ysize=400
	plot,xx,yy,yrange=[min(yy),max(yy)],title=strmid(datfiles(ii),0,strlen(datfiles(ii))-4),BACKGROUND = 255, COLOR = 0 & oplot,xx,yy, psym=1,SYMSIZE=1&oplot,x1,y1,PSYM=2,SYMSIZE=1, COLOR = 0
;plot, xx,yy,XRANGE=[min(xx),max(xx)],CHARSIZE=1, CHARTHICK=1,TICKLEN=1,YRANGE=[MIN(yy),1.1*MAX(yy)],/XSTYLE,BACKGROUND = 255, COLOR = 0&oplot,x1,y1,PSYM=4,SYMSIZE=2, COLOR = 0

	;;plot, igi_work(0,*),xtickname=xtickv2,igi_work(1,*),xGRIDSTYLE=1,xstyle=1,xtitle=xtit,title=field(b),xrange=[round(igi_work(0,0)),round(igi_work(0,nn-1))],yrange=[min(igi_work(1,*)),max(igi_work(1,*))]
	gifname=strmid(datfiles(ii),0,strlen(datfiles(ii))-4)+'.gif'
	out=bytscl(TVRD())
	tvlct,r,g,b,/get
	;write_gif,gifname,out,r,g,b
	;print,xx,yy
endif else begin
	;x=fltarr(n0+1)
	;y=fltarr(n0+1)
	x=dblarr(n0+1)
	y=dblarr(n0+1)
	x(1:n0)=xx(*)
	y(1:n0)=yy(*)
endelse

w=datfiles(ii)
u='d'
D1=5;1
D2=50;N0
L1=0
m=0.0
p=0.0
L=0.0
 PRINT,"MOR6_CO>  6-MORLE ÀÌÏËÈÒÓÄÍÎ-×ÀÑÒÎÒÍÛÉ ÏÎÐÒÐÅÒ È ÓÐÎÂÅÍÜ ÇÍÀ×ÈÌÎÑÒÈ"
 PRINT, " +--Ð å æ è ì û :--------------------------"
 PRINT, " | "
 PRINT, " | N = ",N0
 PRINT, " | * Èìÿ ôàéëà äëÿ ÂÂÎÄà ",w
 PRINT, " | Äëÿ åãî çàïèñè â èìåíè èñïîëüçîâàòü +áóêâó ",u
 PRINT, " | Èíòåðâàëû ïåðèîäîâ (îò N1, äî N2)",d1,d2
 D3=D2-D1+1.0
 Q6=1
 Q4=1
 B6=1
 PRINT, " +-----------------------------------------"
 ;;X=FLTARR(2*N0+1)&Y=fLTARR(N0+1)&
 ;Z=FLTARR(N0+1)&F=FLTARR(2*N0+1)&F1=FLTARR(2*N0+1)&
 ;R=FLTARR(2*N0+1)&LL=FLTARR(50+1)&Q1=FLTARR(2*N0+1)&Q=FLTARR(2*N0+1)&R1=FLTARR(2*N0+1)&T=FLTARR(2*N0+1)
 Z=dblARR(N0+1)&F=DBLARR(2*N0+1)&F1=DBLARR(2*N0+1)&
 R=DBLARR(2*N0+1)&LL=DBLARR(50+1)&Q1=DBLARR(2*N0+1)&Q=DBLARR(2*N0+1)&R1=DBLARR(2*N0+1)&T=DBLARR(2*N0+1)

 OPENR,1, W; FOR INPUT AS #1

 ;FOR I=1, N0 DO BEGIN
 ;readf,1, qq,ww,FORMAT = '(f3,f8)'
 ;z(i)=qq
 ;y(i)=ww
 ;ENDFOR
 ;print,z
 ;print,y
 ;CLOSE,1
z=x
 Q7=Z(1)
; W=U+W
w=strmid(w,0,strlen(w)-4)+u+'.dat'
 close,1
 OPENW,1, W
 CLOSE,1
 goto, LAB510

 A1=0.0
 A2=0.0
 A3=0.0
 A4=0.0
 A5=0.0
 for i=3, N0-2 DO BEGIN
 A1=A1+(I-(N0+1.0)/2.0)
 A2=A2+(I-(N0+1.0)/2.0)*(I-(N0+1)/2.0)
 A3=A3+Y(I)*(I-(N0+1.0)/2.0)
 A4=A4+Y(I)
 A5=A5+Y(I)*Y(I)
ENDFOR
 N0=N0-4.0
 F2=(N0*A3-A1*A4)/(N0*A2-A1*A1)
 G2=(A4-F2*A1)/N0
 S2=SQRT((A5-G2*A4-F2*A3)/(N0-2))
 F3=S2*sqrT(N0/(N0*A2-A1*A1))
 G3=S2*sqrT(A2/(N0*A2-A1*A1))
 ;print,f2,f3,g3
 N0=N0+4.0
 PRINT, "A =",G2,"+-",G3,"     B =",F2/8*60,"+-",F3/8*60,"[1/hour]"
 FOR I=1, N0 DO BEGIN
 F1(I)=Y(I)-G2-F2*(I-(N0+1.0)/2.0)
 ENDFOR
 goto, LAB526


 LAB510:
 A1=0.0
 yy=y
 FOR I=1, N0 DO BEGIN
 A1=A1+Y(I)
 ENDFOR
 U1=A1/N0
 FOR I=1, N0 DO BEGIN
 F1(I)=Y(I)-U1
 ;print, f1(i),y(i),u1,a1
 ENDFOR
 ;print,y
 ;print, a1,u1
 LAB526:
 U1=0.0
 ;;H=FLTARR(D3,N0)
 H=FLTARR(D3+1,N0+1)
 H=dblARR(D3+1,N0+1)

 A1=13.6
 ;REM PRINT "* Periods:"
 FOR I=1, 30 DO BEGIN
 ;REM A1=A1*1.135
 ;REm A2=FIX((A1)/2*2)
 ;REM IF A2/2=A2\2 THEN L(I)=A2+1+(Q4-2)*10 ELSE L(I)=A2+(Q4-2)*10
 ;REM IF i=4 THEN L(I)=11
 ;REM L(I)=FIX(2^(((I+3)/4+1.25)))+1
 ;L(I)=INT(2^((I+7)/4)+.5)+1
 LL(I)=floor(2.0^((I+7.0)/4.0)+0.5)+1.0
 ;print,LL(i)
 ;REM PRINT L(I)-1;
 ENDFOR;NEXT
 M=L1*15.0
 ;M=LL(I)*15.0
 ;print,m

 FOR J=D1+M, D2+M DO BEGIN
 IF LL(J) GT N0 THEN GOTO, LAB1345
 N=LL(J)-1.0
 ;PRINT, FORMAT = '( 12(i5),"-")', J

	 FOR I=-N, N DO BEGIN
 	;D=2.0*2.0*I*I/N/N
 	D=(4*I*I)/N/N
 	R1(I+N)=COS(2.0*6.0*I/N)*EXP(-D/2.0)/SQRT(N)
 	Q1(I+N)=SIN(2.0*6.0*I/N)*EXP(-D/2.0)/SQRT(N)
 	;print, r1(i+n),q1(i+n)
	ENDFOR;NEXT I

	 ;FOR K=N*2\6+1, N0-N*2\6 DO BEGIN
	 kbeg=fix(N*2.0/6.0)+1.0
	 kend=N0-fix(N*2.0/6.0)
	 ;print, kbeg,kend
 	;FOR K=fix(N*2.0/6.0)+1, N0-fix(N*2.0/6.0) DO BEGIN
 	FOR K=kbeg, kend DO BEGIN
 		FOR I=-N, N DO BEGIN
 		IF ((I+K) LT 0) OR ((I+K) GT N0) THEN GOTO, LAB816 ELSE GOTO, LAB821
 		LAB816:
 		F(I+N)=U1
 		R(I+N)=0.0
 		Q(I+N)=0.0
 		goto, LAB839
 		LAB821:
 		F(I+N)=F1(I+K)
 		R(I+N)=R1(I+N)
 		Q(I+N)=Q1(I+N)
 	;	print, f(i+n),r(i+n),q(i+n)
 		LAB839:
 	;;;	print, f(i+n),r(i+n),q(i+n)
	 	ENDFOR;NEXT I
;print, f,r,q
	 	A8=0.0
 		C8=0.0
 		A5=0.0
		C5=0.0

 		FOR I=-N, N DO BEGIN
 		A8=A8+R(I+N)*F(I+N)
 		C8=C8+Q(I+N)*F(I+N)
 		A5=A5+R(I+N)*R(I+N)
 		C5=C5+Q(I+N)*Q(I+N)
 		ENDFOR;NEXT I
 		;print,a8,c8,a5,c5
	 B8=A8/A5
 	V8=C8/C5
 	;REM B8=A8/SQR(A5*A5+C5*C5)
 	;REM V8=C8/SQR(A5*A5+C5*C5)
 	B8=A8
 	V8=C8
 	;IF ABS(B8)LT 0.001 THEN T4=0 ELSE T4=CSNG(INT(B8*1000000+.5)/1000000)
 	;IF ABS(V8) LT 0.001 THEN T5=0 ELSE T5=CSNG(INT(V8*1000000+.5)/1000000)
 	;IF ABS(B8)LT 0.001 THEN T4=0.0 ELSE T4=float(fix(B8*1000000.0+0.5)/1000000.0)
 	;IF ABS(V8) LT 0.001 THEN T5=0.0 ELSE T5=float(fix(V8*1000000.0+0.5)/1000000.0)

 	IF ABS(B8) LT 0.001 THEN T4=0.0 ELSE T4=float(long(B8*1000000+0.5)/1000000)
 	IF ABS(V8) LT 0.001 THEN T5=0.0 ELSE T5=float(long(V8*1000000+0.5)/1000000)
 	;;;print, b8, v8, t4,t5
 	N6=J-M-D1+1
 	;H(N6,K)=CSNG(INT(SQRT(T4*T4+T5*T5)*2*sqrT(2)/sqrT(3.14159*N)*10000+.5)/10000)
 	H(N6,K)=float(long(SQRT(T4*T4+T5*T5)*2*sqrT(2)/sqrT(3.14159*N)*10000.0+0.5)/10000.0)
 	;REM H(N6,K)=T4
 	ENDFOR;NEXT K

 ;PRINT, "-";
 ENDFOR;NEXT J

 LAB1345:
 PRINT, "* Central Periods:"
 FOR I=1, 30 DO BEGIN
 ;REM A1=A1*1.135
 ;REM A2=FIX(A1)
 ;REM IF A2/2=A2\2 THEN L(I)=A2+1+(Q4-2)*10 ELSE L(I)=A2+(Q4-2)*10
 ;REM IF I=4 THEN L(I)=11
 ;REM L(I)=FIX(2^(((I+3)/4+1.25)))+1
 ;REM L(I)=INT(2^((I+7)/4)+.5)+1
 ;REM PRINT L(I)-1;
 X(I)=(LL(I)-1.0)*2.0*3.14159/12.0
 ;PRINT, INT(X(I)*10+0.5)/10;
 ;PRINT, fix(X(I)*10.0+0.5)/10.0; printprint
 ENDFOR;NEXT
;PRINT, fix(X(1:30)*10.0+0.5)/10.0
 ;PRINT
 ;PRINT
 PRINT, "____________________________"
 OPENW,1, W;$ FOR OUTPUT AS #1
 CLOSE,1

 FOR J=1, N0 DO BEGIN
	 FOR I=D1, D2 DO BEGIN
 	 Y(I)=H(I-D1+1.0,J)
	 ENDFOR;NEXT I
	 ;print,J

 	LAB1512:
 	FOR I=D1+1, D2-1 DO BEGIN
 	;print,I
 		IF Y(I+1) EQ 0 THEN GOTO, LAB1670
 		IF (Y(I) gt Y(I-1)) AND (Y(I) GT Y(I+1)) THEN GOTO, LAB1520 ELSE GOTO, LAB1640
		LAB1520:
 		m1=X(I-1)*X(I-1)*X(I)-X(I-1)*X(I)*X(I)
 		m2=X(I-1)*X(I-1)-X(I)*X(I)
 		m3=Y(I)*X(I-1)*X(I-1)-Y(I-1)*X(I)*X(I)
 		m4=X(I)*X(I)*X(I+1)-X(I)*X(I+1)*X(I+1)
 		m5=X(I)*X(I)-X(I+1)*X(I+1)
 		m6=Y(I+1)*X(I)*X(I)-Y(I)*X(I+1)*X(I+1)
 		C=(m1*m6-m3*m4)/(m1*m5-m2*m4)
 		B=(m3-C*m2)/m1
 		A=(Y(I)-C-B*X(I))/X(I)/X(I)
 		S=-B/2.0/A
 	 	M=A*S*S+B*S+C
 		IF (S GT 0.2) AND (S LT 400.0) THEN GOTO, LAB1532 ELSE GOTO, LAB1636
	 	LAB1532:
	 	;print,'1532'
 		;IF (j/((L(I)-1)\24+1)) EQ (j\((L(I)-1)\24+1)) THEN GOTO, LAB1533 ELSE 	GOTO, LAB1636
 		IF (j/(FIX((LL(I)-1.0)/24.0)+1)) EQ FIX(j/(FIX((LL(I)-1.0)/24.0)+1.0))$
 		THEN GOTO, LAB1533 ELSE GOTO, LAB1636
		LAB1533:;----------------------
		;print, '1533'
 		N=S*12.0/2.0/3.14159
 		FOR L=-N, N DO BEGIN
 		;print,L
 		D=2.0*2.0*L*L/N/N
 	    R1(L+N)=COS(2.0*6.0*L/N)*EXP(-D/2.0)/SQRT(N)
 		Q1(L+N)=SIN(2.0*6.0*L/N)*EXP(-D/2.0)/SQRT(N)
 		ENDFOR;NEXT L

		 K=j
	 	LAB1542:;------------------------
	 	;print,'1542'
	  	FOR L=-N, N DO BEGIN
	 	;print,L
	 	IF ((L+K) LT 0)OR((L+K) GT N0) THEN GOTO, LAB1546 ELSE GOTO, LAB1554
		LAB1546:
 		F(L+N)=U1
 		R(L+N)=0.0
 		Q(L+N)=0.0
	 	goto, LAB1560
 		LAB1554:
 		F(L+N)=F1(L+K)
 		R(L+N)=R1(L+N)
 		Q(L+N)=Q1(L+N)
 		LAB1560:
 		ENDFOR;NEXT L

 		A8=0.0
 		C8=0.0
 		A5=0.0
 		C5=0.0

 		FOR L=-N, N DO BEGIN
 		A8=A8+R(L+N)*F(L+N)
 		C8=C8+Q(L+N)*F(L+N)
 		A5=A5+R(L+N)*R(L+N)
 		C5=C5+Q(L+N)*Q(L+N)
 		ENDFOR;NEXT L

 		B8=A8/A5
 		V8=C8/C5
	 	;REM B8=A8/SQRT(A5*A5+C5*C5)
		 ;REM V8=C8/SQRT(A5*A5+C5*C5)
 		B8=A8
 		V8=C8
 		;;IF ABS(B8) LT 0.001 THEN T4=0 ELSE T4=CSNG(INT(B8*1000000+.5)/1000000)
 		;;IF ABS(V8) LT 0.001 THEN T5=0 ELSE T5=CSNG(INT(V8*1000000+.5)/1000000)
 		;IF ABS(B8) LT 0.001 THEN T4=0.0 ELSE T4=float(long(B8*1000000.0+0.5)/1000000.0)
 		;IF ABS(V8) LT 0.001 THEN T5=0.0 ELSE T5=float(long(V8*1000000.0+0.5)/1000000.0)
		IF ABS(B8) LT 0.001 THEN T4=0.0 ELSE T4=double(long(B8*1000000.0+0.5)/1000000.0)
 		IF ABS(V8) LT 0.001 THEN T5=0.0 ELSE T5=double(long(V8*1000000.0+0.5)/1000000.0)

 		FOR L=-N, N DO BEGIN
 		T(L+N)=T4*R(L+N)+T5*Q(L+N)
 		ENDFOR;NEXT L

 		A8=0.0
 		C8=0.0
 		A5=0.0
 		C5=0.0
 		A6=0.0

	 	FOR L=-N, N DO BEGIN
 		A8=A8+T(L+N)
 		C8=C8+T(L+N)*T(L+N)
	 	A5=A5+T(L+N)*F(L+N)
	 	C5=C5+F(L+N)
 		A6=A6+F(L+N)*F(L+N)
 		ENDFOR;NEXT L
close,2

 		P=((2.0*N+1.0)*A5-A8*C5)^2/((2.0*N+1.0)*C8-A8*A8)/((2.0*n+1.0)*A6-C5*C5)
 		P=1.0-(1.0-P)^((2.0*N+1.0-2.0)/2.0-1.0)
 		OPENU,2, W,/APPEND;$ FOR APPEND AS #2
 		PRINTF, 2, Z(1)+(J-1)*(Z(2)-Z(1)), 100.0/S/(Z(2)-Z(1)), M, P
 		;PRINT, 'file: ',Z(1)+(J-1)*(Z(2)-Z(1)), 100.0/S/(Z(2)-Z(1)), M, P
 		CLOSE,2
 		LAB1636:
 		;print,'1636'
 		GOTO, LAB1670
 		LAB1640:
 		IF (Y(I) EQ Y(I-1))AND(Y(I) EQ Y(I+1)) THEN GOTO, LAB1542 ELSE GOTO, LAB1670
 		IF (Y(I) EQ 0)AND(Y(I-1) EQ 0)AND(Y(I+1) EQ 0) THEN GOTO, LAB1670 ELSE GOTO, LAB1647
	 	LAB1647:
 		S=X(I)
 		M=Y(I)
 		OPENU,3, W,/APPEND;$ FOR APPEND AS #3
 		PRINTF,3, J, 100.0/S, M
 		;PRINT, 'file2:', J, 100.0/S, M
 		;PRINT, J, 100.0/S, M
 		CLOSE,3
 		LAB1670:
 		;print,'1670'
	 	K=K
 	ENDFOR;NEXT I
 ENDFOR;NEXT J
 ;PRINT, 'file: ',Z(1)+(J-1)*(Z(2)-Z(1)), 100.0/S/(Z(2)-Z(1)), M, P

 ;PRINT, "---> END"

out=read_ascii(w,count=N1)
;col=fltarr(4,n1)
col=dblarr(4,n1)
window,0,xsize=600,ysize=1000
xxx=out.field1(0,*);)/max(out.field1(0,*))
yyy=100.0/out.field1(1,*);)*400)/max(100.0/out.field1(1,*))
plot,xxx,yyy,xstyle=1,title=strmid(datfiles(ii),0,strlen(datfiles(ii))-4),xrange=[round(min(xxx)),round(max(xxx))],ystyle=1,yrange=[min(yyy),max(yyy)],psym=6,BACKGROUND = 255, COLOR = 0;,SYMSIZE=out.field1(3,*)/10.0
;out.field01(0,*) out.field01(1,*)
xxx1=(out.field1(0,*)-min(out.field1(0,*)))
xxx1=xxx1/max(xxx1)
xxx1=xxx1*510+70
yyy11=((100.0/out.field1(1,*)-min(100.0/out.field1(1,*))))
yyy11=yyy11/max(yyy11)
yyy1=yyy11*320+35
zzz=out.field1(2,*)*8/max(out.field1(2,*))


for kk=0, n0 do begin
XYOUTS, xxx1(kk), yyy1(kk), 'o', CHARSIZE=zzz(kk),ALIGNMENT=0.5, COLOR = 0, /DEVICE
endfor

;plot,xxx,yyy,yrange=[min(yyy),max(yyy)],psym=6;,SYMSIZE=out.field1(3,*)/10.0 ;& oplot,xx,yy, psym=2,SYMSIZE=2&oplot,x1,y1,PSYM=4,SYMSIZE=2
;plot, igi_work(0,*),xtickname=xtickv2,igi_work(1,*),xGRIDSTYLE=1,xstyle=1,xtitle=xtit,title=field(b),xrange=[round(igi_work(0,0)),round(igi_work(0,nn-1))],yrange=[min(igi_work(1,*)),max(igi_work(1,*))]
gifname=strmid(w,0,strlen(w)-4)+'w.gif'
out=bytscl(TVRD())
tvlct,r,g,b,/get
write_gif,gifname,out,r,g,b

kon1:
endfor
kon:
CLOSE,1
 CLOSE,2
 CLOSE,3
 PRINT, "---> END"
end