pro rectangle,xx,yy,sirka,vyska ;nakresli obdelnik ;device,set_graphics=6 px=[xx,xx+sirka,xx+sirka,xx,xx] py=[yy,yy,yy+vyska,yy+vyska,yy] plots,px,py,/dev,color=255;lines=0,color=255 ;empty end pro rectangle1,xx,yy,sirka,vyska ;device,set_graphics=6 px=[xx,xx+sirka,xx+sirka,xx,xx] py=[yy,yy,yy+vyska,yy+vyska,yy] plots,px,py,/dev,color=0;lines=0,color=255 ;empty end function integruj,obr,x,y,w,h ;spocita integral v oblasti x:x+w,y:y+h sub=float(obr(x:x+w,y:y+h)) & s=size(sub) iks=findgen(s(1)) & ips=findgen(s(2)) pole=float(findgen(s(2))) for i=0,s(2)-1 do pole(i)=int_tabulated(iks,sub(*,i)) return,int_tabulated(ips,pole) end pro hist2_event,event COMMON uloz4,histdraw2,corrdraw2,corrbase2 widget_control,event.id,get_uvalue=uvalue case uvalue of 'SAVE':begin widget_control,histdraw2,get_value=okno2 wset,okno2 & write_jpeg,'hist1.jpg',tvrd() if widget_info(corrbase2,/valid_id) eq 1 then begin widget_control,corrdraw2,get_value=okno1 wset,okno1 & write_jpeg,'corr1.jpg',tvrd() widget_control,corrbase2,/destroy endif widget_control,event.top,/destroy end 'CANCEL':begin widget_control,event.top,/destroy if widget_info(corrbase2,/valid_id) eq 1 then widget_control,corrbase2,/destroy end endcase end pro cross_event,event COMMON uloz3,crossdraw2 widget_control,event.id,get_uvalue=uvalue case uvalue of 'OK':begin widget_control,crossdraw2,get_value=okno & wset,okno write_jpeg,'cross.jpg',tvrd() widget_control,event.top,/destroy end 'CANCEL':widget_control,event.top,/destroy endcase end pro auto_event,event COMMON uloz2,autodraw2 widget_control,event.id,get_uvalue=uvalue case uvalue of 'OK':begin widget_control,autodraw2,get_value=okno & wset,okno write_jpeg,'auto.jpg',tvrd() widget_control,event.top,/destroy end 'CANCEL':widget_control,event.top,/destroy endcase end pro light_event,event COMMON uloz6,lightdraw2 widget_control,event.id,get_uvalue=uvalue case uvalue of 'OK':begin widget_control,lightdraw2,get_value=okno & wset,okno write_jpeg,'light.jpg',tvrd() widget_control,event.top,/destroy end 'CANCEL':widget_control,event.top,/destroy endcase end pro hist_event,event COMMON uloz,corrdraw,histdraw widget_control,event.id,get_uvalue=uvalue case uvalue of 'OK':begin widget_control,corrdraw,get_value=okno1 widget_control,histdraw,get_value=okno2 wset,okno1 & write_jpeg,'corr.jpg',tvrd() wset,okno2 & write_jpeg,'hist.jpg',tvrd() widget_control,event.top,/destroy end 'CANCEL':widget_control,event.top,/destroy endcase end pro mnoho_event,event COMMON uloz5,mnohofield,mnohodraw,mnohodraw2 widget_control,event.id,get_uvalue=uvalue case uvalue of 'OK':begin widget_control,mnohodraw,get_value=okno1 widget_control,mnohodraw2,get_value=okno2 wset,okno1 & write_jpeg,'reg.jpg',tvrd() wset,okno2 & write_jpeg,'count.jpg',tvrd() widget_control,event.top,/destroy end 'CANCEL':widget_control,event.top,/destroy endcase end ;--------zde se pocitaji integraly pro vsechny oblasti site pro kazdy obrazek pro integrals_event,event COMMON integral,poleintegralu,xpocatek,ypocatek,xpocet,ypocet,krok,prumer COMMON obrazky,soubory,pocet COMMON oblast,inttext1,inttext2,inttext3 COMMON baze,base,innerbase widget_control,event.id,get_uvalue=uvalue case uvalue of 'OK':begin ;-----------nacteni zadani oblasti---------------------- widget_control,inttext1,get_value=text1 & text1=text1(0) widget_control,inttext2,get_value=text2 & text2=text2(0) widget_control,inttext3,get_value=text3 & text3=text3(0) pozice=strpos(text1,',') & delka=strlen(text1) xpocatek=fix(strmid(text1,0,pozice)) & ypocatek=fix(strmid(text1,pozice+1,delka-pozice)) pozice=strpos(text2,',') & delka=strlen(text2) xpocet=fix(strmid(text2,0,pozice)) & ypocet=fix(strmid(text2,pozice+1,delka-pozice)) krok=fix(text3) widget_control,event.top,/destroy ;----------vytvoreni pole integralu-------------------- obrazek=readfits(soubory(pocet/2)) & s=size(obrazek) widget_control,innerbase,/destroy innerbase=widget_base(base,/column) d=widget_draw(innerbase,xsize=s(1),ysize=s(2)) innerlabel=widget_label(innerbase,xsize=900,/align_left) poleintegralu=fltarr(xpocet,ypocet,pocet) for o=0,pocet-1 do begin ;pocet-1 obrazek=readfits(soubory(o)) widget_control,innerlabel,set_value=soubory(o) for i=0,xpocet-1 do begin for j=0,ypocet-1 do begin poleintegralu(i,j,o)=integruj(obrazek,xpocatek+i*krok,ypocatek+j*krok,krok,krok) endfor endfor endfor ;pro dany obrazek odecteme od integralu jejich stredni hodnotu prumer=fltarr(pocet) for i=0,pocet-1 do begin ;pocet-1 prumer(i)=total(poleintegralu(*,*,i))/float(xpocet*ypocet) poleintegralu(*,*,i)=poleintegralu(*,*,i)-prumer(i) endfor ;ulozeni integralu, strednich hodnot a parametru site a vykresleni site save,prumer,filename='average.dat' save,poleintegralu,filename='variable.dat' openw,unit,'region.dat',/get_lun printf,unit,xpocatek,ypocatek,xpocet,ypocet,krok ;ulozeni parametru oblasti free_lun,unit tvscl,obrazek for i=0,xpocet-1 do begin for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok endfor widget_control,innerlabel,set_value='The integrals have been counted' end 'CANCEL':widget_control,event.top,/destroy endcase end pro korel_event,event COMMON uloz4,histdraw2,corrdraw2,corrbase2 COMMON baze,base,innerbase COMMON obrazky,soubory,pocet COMMON oblast,inttext1,inttext2,inttext3 COMMON integral,poleintegralu,xpocatek,ypocatek,xpocet,ypocet,krok,prumer COMMON start,omenu COMMON uloz,corrdraw,histdraw COMMON uloz2,autodraw2 COMMON uloz3,crossdraw2 COMMON inner1,autofield COMMON inner2,crossfield1,crossfield2 COMMON inner3,distfield COMMON inner4,corectfield1,corectfield2,corectfield3 COMMON inner5,prepinac,prepinac1 COMMON uloz5,mnohofield,mnohodraw,mnohodraw2 COMMON inner6,lightfield COMMON uloz6,lightdraw2 widget_control,event.id,get_uvalue=uvalue case uvalue of 'LOAD':begin ;zde se nahrava adresar se serii obrazku adresar=dialog_pickfile(/directory) soubory=findfile(adresar+'*.fts',count=pocet) ;pokud jsou spocteny integraly, nahrajem je r=findfile('variable.dat') ;integraly if r(0) ne '' then restore,'variable.dat' reg=findfile('region.dat') ;parametry oblasti if reg(0) ne '' then begin openr,unit,'region.dat',/get_lun readf,unit,xpocatek,ypocatek,xpocet,ypocet,krok free_lun,unit endif pr=findfile('average.dat') ;stredni hodnoty if pr(0) ne '' then restore,'average.dat' widget_control,omenu,sensitive=1 end 'QUIT':widget_control,event.top,/destroy ;----------------------------------------------------------- 'INTEGRALS':begin intbase=widget_base(/column,title='region for integration') intbase1=widget_base(intbase,/column,/frame) inttext1=cw_field(intbase1,title='Start point (x,y)') inttext2=cw_field(intbase1,title='Number of regions (n,m)') inttext3=cw_field(intbase1,title='Step (k)') intbase2=widget_base(intbase,/row) okbut=widget_button(intbase2,value='OK',uvalue='OK') canbut=widget_button(intbase2,value='CANCEL',uvalue='CANCEL') widget_control,intbase,/realize xmanager,'integrals',intbase end 'KOEFICIENTY':begin ;zde se pocitaji korelacni koeficienty pro vsechny dvojice oblasti obrazek=readfits(soubory(pocet/2)) & s=size(obrazek) widget_control,innerbase,/destroy innerbase=widget_base(base,/column) d=widget_draw(innerbase,xsize=s(1),ysize=s(2)) tvscl,obrazek ;------------------ukladani korelacnich koeficientu----------- openw,unit,'koeficienty.dat',/get_lun for i=0,ypocet-1 do begin ;kombinace oblasti for j=0,xpocet-1 do begin ;oblast [l,k] a [j,i] for k=0,ypocet-1 do begin for l=0,xpocet-1 do begin if (k gt i) or ((k eq i) and (l gt j)) then begin x1=xpocatek+l*krok & y1=ypocatek+k*krok x2=xpocatek+j*krok & y2=ypocatek+i*krok lag=0 & x=poleintegralu(l,k,*) & y=poleintegralu(j,i,*) result=c_correlate(x,y,lag) printf,unit,result(0),x2,y2,x1,y1 rectangle,x2,y2,krok,krok endif endfor endfor endfor endfor free_lun,unit corlabel=widget_label(innerbase,value='The correlation coefficients has been saved') end 'HISTOGRAM':begin ;zde se pocita histogram vsech korelacnich koeficientu widget_control,innerbase,/destroy innerbase=widget_base(base,/column) obrazek=readfits(soubory(pocet/2)) s=size(obrazek) corrdraw=widget_draw(innerbase,xsize=s(1),ysize=s(2)) tvscl,obrazek ;----------nacteni koeficientu ze souboru-------------- n=long(xpocet)*long(ypocet) ;nezapomenout na longint! pocetkombinaci=n*(n-1)/2 openr,unit,'koeficienty.dat',/get_lun data=fltarr(5,pocetkombinaci) readf,unit,data free_lun,unit ;---------rozdeleni a jeho parametry-------- hist=histogram(data(0,*),min=-1,max=1,bins=0.01) hodnotyx=findgen(n_elements(hist))*0.01-1 gauss=gaussfit(hodnotyx,hist,koef,nterms=3) mi=koef(1) & sigma=koef(2) horni=mi+3.0*sigma & dolni=mi-3.0*sigma ;---------vykresleni oblasti s vysokymi korelacemi------- indexy=where((data(0,*) lt dolni) or (data(0,*) gt horni),count) korelace=data(*,indexy) for i=0,count-1 do begin rectangle,korelace(1,i),korelace(2,i),krok,krok ;x,y,s,v rectangle,korelace(3,i),korelace(4,i),krok,krok endfor ;--------zobrazeni histogramu-------------- histbase=widget_base(/column) histdraw=widget_draw(histbase,xsize=500,ysize=400) histbase1=widget_base(histbase,/row) okbut=widget_button(histbase1,value='Save images',uvalue='OK') canbut=widget_button(histbase1,value='Cancel',uvalue='CANCEL') widget_control,histbase,/realize xmanager,'hist',histbase widget_control,histdraw,get_value=okno2 & wset,okno2 plot,hodnotyx,hist,xrange=[-1,1],xstyle=1,$ title='DISTRIBUTION OF CORRELATION COEFFICIENTS',$ xtitle='correlation coefficient',ytitle='density' oplot,hodnotyx,gauss,linestyle=1 xyouts,100,303,'...........',/device xyouts,150,300,'gauss fit',/device xyouts,100,280,'mi='+strcompress(string(mi),/remove_all),/device xyouts,100,265,'sigma='+strcompress(string(sigma),/remove_all),/device end 'MNOHO':begin widget_control,innerbase,/destroy innerbase=widget_base(base,/row) obrazek=readfits(soubory(pocet/2)) s=size(obrazek) mnohodraw=widget_draw(innerbase,xsize=s(1),ysize=s(2)) mnohobase=widget_base(innerbase,/column,xsize=250) mnohofield=cw_field(mnohobase,title='n ') mnohobutton=widget_button(mnohobase,value='SHOW REGIONS',uvalue='MNOHO2') end 'MNOHO2':begin ;zde se zjistuje s kolika dalsimi oblastmi oblasti koreluji widget_control,mnohofield,get_value=kolik & kolik=fix(kolik(0)) widget_control,mnohodraw,get_value=okno & wset,okno erase ;------nacteni korelacnich koeficientu ze souboru------- n=long(xpocet)*long(ypocet) pocetkombinaci=n*(n-1)/2 openr,unit,'koeficienty.dat',/get_lun data=fltarr(5,pocetkombinaci) readf,unit,data ;zde jsou vsechny korelacni koeficienty free_lun,unit ;-----spocteni histogramu a zjisteni mi a sigma--------- hist=histogram(data(0,*),min=-1,max=1,bins=0.01) hodnotyx=findgen(n_elements(hist))*0.01-1 gauss=gaussfit(hodnotyx,hist,koef,nterms=3) mi=koef(1) & sigma=koef(2) horni=mi+3.0*sigma & dolni=mi-3.0*sigma ;------z pole data vytrhneme vysoke korelace----------- indexy=where((data(0,*) lt dolni) or (data(0,*) gt horni),count) korelace=data(*,indexy) ;pole s vysokymi korelacemi ;----vytvorime pole oblasti, ve kterem budou pocty oblasti, s kterymi ;----prislusna oblast koreluje--------------------- oblasti=intarr(xpocet,ypocet) for i=0,count-1 do begin x1=korelace(1,i) & y1=korelace(2,i) x2=korelace(3,i) & y2=korelace(4,i) s1=(x1-xpocatek)/krok & s2=(y1-ypocatek)/krok oblasti(s1,s2)=oblasti(s1,s2)+1 s1=(x2-xpocatek)/krok & s2=(y2-ypocatek)/krok oblasti(s1,s2)=oblasti(s1,s2)+1 endfor ;----vykresleni oblasti, ktere koreluji s vice nez n oblastmi obrazek=readfits(soubory(pocet/2)) tvscl,obrazek for i=0,xpocet-1 do begin for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok endfor for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5 for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5 indexy1=where(oblasti gt kolik,count1) if count1 ne 0 then begin print,'souradnice oblasti:' for i=0,count1-1 do begin s1=fix(indexy1(i) mod xpocet) s2=fix(indexy1(i)/xpocet) rectangle1,s1*krok+xpocatek,s2*krok+ypocatek,krok,krok ;print,s1+1,s2+1 endfor endif ;---------pocty oblasti a jejich contour------------------ mnohobase2=widget_base(/column) mnohodraw2=widget_draw(mnohobase2,xsize=400,ysize=600) mnohobase21=widget_base(mnohobase2,/row) mnohobut1=widget_button(mnohobase21,value='SAVE',uvalue='OK') mnohobut2=widget_button(mnohobase21,value='CANCEL',uvalue='CANCEL') widget_control,mnohobase2,/realize xmanager,'mnoho',mnohobase2 widget_control,mnohodraw2,get_value=okno & wset,okno obr=obrazek(xpocatek:xpocatek+xpocet*krok,ypocatek:ypocatek+ypocet*krok) obr2=obr ;-----------------pocty oblasti------------- mez=50 konst=fix(float(mez*max(oblasti))/float(255-mez)) for m=0,xpocet-1 do begin for n=0,ypocet-1 do begin xx=m*krok & yy=n*krok jas=fix(255*float(oblasti(m,n)+konst)/float(max(oblasti)+konst)) if oblasti(m,n) eq 0 then jas=0 obr2(xx:xx+krok,yy:yy+krok)=jas ;zde oblasti priradime jas podle poctu endfor ;oblasti s kterymi koreluje endfor tv,congrid(obr2,370,270),15,315 for m=0,xpocet-1 do begin xkrok=370.0/xpocet & ykrok=270.0/ypocet for n=0,ypocet-1 do rectangle,m*xkrok+15,n*ykrok+315,xkrok,ykrok endfor xyouts,200,590,'COUNT OF CORRELATIONS',charsize=1,alignment=0.5,/device ;----------------contour------------------ tvscl,congrid(obr,370,270),15,15 contour,oblasti,charsize=0.7,title='CONTOUR',pos=[15,15,385,285],/noerase,xmargin=[15,15],$ ymargin=[15,15],xrange=[-0.5,xpocet-1.5],yrange=[0.5,ypocet-1.5],xstyle=1,ystyle=1,$ xticks=xpocet,yticks=ypocet,xtickn=replicate(' ',xpocet+1),ytickn=replicate(' ',ypocet+1),/device end 'AUTO':begin ;----------zobrazeni obrazku--------------- widget_control,innerbase,/destroy innerbase=widget_base(base,/row) obrazek=readfits(soubory(pocet/2)) s=size(obrazek) autodraw=widget_draw(innerbase,xsize=s(1),ysize=s(2)) tvscl,obrazek ;-------------vykresleni site----------- for i=0,xpocet-1 do begin for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok endfor for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5 for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5 autobase=widget_base(innerbase,/column,xsize=250) autofield=cw_field(autobase,title='coordinates of region (x,y)') autobut=widget_button(autobase,value='CORRELATE',uvalue='AUTOCORR') end 'AUTOCORR':begin ;zde se pocita autokorelace vybrane oblasti a jeji FFT widget_control,autofield,get_value=sourad & sourad=sourad(0) pozice=strpos(sourad,',') & delka=strlen(sourad) xsourad=fix(strmid(sourad,0,pozice)) & ysourad=fix(strmid(sourad,pozice+1,delka-pozice)) ;----graf s korelaci-------- autobase2=widget_base(/column) autodraw2=widget_draw(autobase2,xsize=400,ysize=600) autobase3=widget_base(autobase2,/row) autobut2=widget_button(autobase3,value='SAVE',uvalue='OK') autocan=widget_button(autobase3,value='CANCEL',uvalue='CANCEL') widget_control,autobase2,/realize xmanager,'auto',autobase2 widget_control,autodraw2,get_value=okno & wset,okno rada=poleintegralu(xsourad-1,ysourad-1,*) lag=findgen(2*pocet-1)-pocet+1 res=a_correlate(rada,lag) koef=res(pocet-1) & koef=strcompress(string(koef),/remove_all) stary=!p.multi & !p.multi=[0,1,2] plot,lag,res,title='AUTO CORRELATION c0='+koef,xtitle='lag',ytitle='correlation coefficient',$ yrange=[-1,1],ystyle=1 plot,/ylog,abs(fft(res,-1)),title='THE LOG OF THE POWER SPECTRUM' !p.multi=stary end 'CROSS':begin widget_control,innerbase,/destroy innerbase=widget_base(base,/row) obrazek=readfits(soubory(pocet/2)) s=size(obrazek) crossdraw=widget_draw(innerbase,xsize=s(1),ysize=s(2)) tvscl,obrazek ;------------zobrazeni site--------------- for i=0,xpocet-1 do begin for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok endfor for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5 for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5 crossbase=widget_base(innerbase,/column,xsize=250) crossfield1=cw_field(crossbase,title='coordinates of 1st region (x,y)') crossfield2=cw_field(crossbase,title='coordinates of 2nd region (x,y)') crossbut=widget_button(crossbase,value='CORRELATE',uvalue='CROSSCORR') end 'CROSSCORR':begin ;zde se pocita korelace dvou vybranych oblasti a jeji FFT widget_control,crossfield1,get_value=sourad & sourad=sourad(0) pozice=strpos(sourad,',') & delka=strlen(sourad) xsourad1=fix(strmid(sourad,0,pozice)) & ysourad1=fix(strmid(sourad,pozice+1,delka-pozice)) widget_control,crossfield2,get_value=sourad & sourad=sourad(0) pozice=strpos(sourad,',') & delka=strlen(sourad) xsourad2=fix(strmid(sourad,0,pozice)) & ysourad2=fix(strmid(sourad,pozice+1,delka-pozice)) ;----graf s korelaci-------- crossbase2=widget_base(/column) crossdraw2=widget_draw(crossbase2,xsize=400,ysize=600) crossbase3=widget_base(crossbase2,/row) crossbut2=widget_button(crossbase3,value='SAVE',uvalue='OK') crosscan=widget_button(crossbase3,value='CANCEL',uvalue='CANCEL') widget_control,crossbase2,/realize xmanager,'cross',crossbase2 widget_control,crossdraw2,get_value=okno & wset,okno rada1=poleintegralu(xsourad1-1,ysourad1-1,*) rada2=poleintegralu(xsourad2-1,ysourad2-1,*) lag=findgen(2*pocet-1)-pocet+1 res=c_correlate(rada1,rada2,lag) koef=res(pocet-1) & koef=strcompress(string(koef),/remove_all) stary=!p.multi & !p.multi=[0,1,2] plot,lag,res,title='CROSS CORRELATION c0='+koef,xtitle='lag',ytitle='correlation coefficient',$ yrange=[-1,1],ystyle=1 plot,/ylog,abs(fft(res,-1)),title='THE LOG OF THE POWER SPECTRUM' !p.multi=stary end 'LIGHT':begin widget_control,innerbase,/destroy innerbase=widget_base(base,/row) obrazek=readfits(soubory(pocet/2)) s=size(obrazek) lightdraw=widget_draw(innerbase,xsize=s(1),ysize=s(2)) tvscl,obrazek ;-----------zobrazeni site--------------- for i=0,xpocet-1 do begin for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok endfor for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5 for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5 lightbase=widget_base(innerbase,/column,xsize=250) lightfield=cw_field(lightbase,title='coordinates of region (x,y)') lightbut=widget_button(lightbase,value='SHOW LIGHT CURVE',uvalue='LIGHT2') end 'LIGHT2':begin ;zde se zobrazi svetelna krivka vybrane oblasti a jeji FFT widget_control,lightfield,get_value=sourad & sourad=sourad(0) pozice=strpos(sourad,',') & delka=strlen(sourad) xsourad=fix(strmid(sourad,0,pozice)) & ysourad=fix(strmid(sourad,pozice+1,delka-pozice)) ;----graf s se svetelnou krivkou a FFT-------- lightbase2=widget_base(/column) lightdraw2=widget_draw(lightbase2,xsize=700,ysize=600) lightbase3=widget_base(lightbase2,/row) lightbut2=widget_button(lightbase3,value='SAVE',uvalue='OK') lightcan=widget_button(lightbase3,value='CANCEL',uvalue='CANCEL') widget_control,lightbase2,/realize xmanager,'light',lightbase2 widget_control,lightdraw2,get_value=okno & wset,okno rada=poleintegralu(xsourad-1,ysourad-1,*) ;rada po odecteni str. hodnoty radapuv=rada+prumer ;rada puvodni - pred odectenim str. hodnoty integralu stary=!p.multi & !p.multi=[0,2,2] plot,rada,title='LIGHT CURVE',charsize=1,xtitle='after substracting the average value',$ yrange=[min(rada)-500,max(rada)+500],ystyle=1 plot,/ylog,abs(fft(rada,-1)),title='THE LOG OF THE POWER SPECTRUM',charsize=1 plot,radapuv,title='LIGHT CURVE',charsize=1,xtitle='before substracting the average value',$ yrange=[min(radapuv)-500,max(radapuv)+500],ystyle=1 plot,/ylog,abs(fft(radapuv,-1)),title='THE LOG OF THE POWER SPECTRUM',charsize=1 !p.multi=stary end 'DIST':begin widget_control,innerbase,/destroy innerbase=widget_base(base,/row) obrazek=readfits(soubory(pocet/2)) s=size(obrazek) crossdraw=widget_draw(innerbase,xsize=s(1),ysize=s(2)) tvscl,obrazek ;-----------zobrazeni site--------------- for i=0,xpocet-1 do begin for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok endfor for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5 for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5 distbase=widget_base(innerbase,/column,xsize=250) distfield=cw_field(distbase,title='coordinates of region (x,y)') togbase1=widget_base(distbase,/column,/exclusive) prep1=widget_button(togbase1,value='histogram from all the coefficients',uvalue='togg1',/no_release) prep2=widget_button(togbase1,value='histogram only for chosen region',uvalue='togg2',/no_release) distbut=widget_button(distbase,value='SHOW CORRELATION',uvalue='DISTCORR') widget_control,togbase1,/realize widget_control,prep1,/set_button prepinac1=1 end 'togg1':prepinac1=1 'togg2':prepinac1=2 'DISTCORR':begin ;zde se zjistuje s kterymi oblastmi zadana oblast koreluje obrazek=readfits(soubory(pocet/2)) s=size(obrazek) corrbase2=widget_base(/column) corrdraw2=widget_draw(corrbase2,xsize=s(1),ysize=s(2)) widget_control,corrbase2,/realize histbase2=widget_base(/column) histdraw2=widget_draw(histbase2,xsize=500,ysize=400) histbase21=widget_base(histbase2,/row) histbut2=widget_button(histbase21,value='SAVE IMAGES',uvalue='SAVE') histcan2=widget_button(histbase21,value='CANCEL',uvalue='CANCEL') widget_control,histbase2,/realize xmanager,'hist2',histbase2 ;----------zadane hodnoty----------------- widget_control,distfield,get_value=sourad & sourad=sourad(0) pozice=strpos(sourad,',') & delka=strlen(sourad) xsourad=fix(strmid(sourad,0,pozice)) & ysourad=fix(strmid(sourad,pozice+1,delka-pozice)) x=(xsourad-1)*krok+xpocatek & y=(ysourad-1)*krok+ypocatek ;----------nacteni korelacnich koeficientu----------------- n=long(xpocet)*long(ypocet) ;nezapomenout na longint!! pocetkombinaci=n*(n-1)/2 openr,unit,'koeficienty.dat',/get_lun data=fltarr(5,pocetkombinaci) readf,unit,data free_lun,unit ;---------vytrhnuti koeficientu tykajicich se zadane oblasti---------- indexy=where(((data(1,*) eq x) and (data(2,*) eq y)) or $ ((data(3,*) eq x) and (data(4,*) eq y)),count) if prepinac1 eq 1 then korelaceh=data else korelaceh=data(*,indexy) korelace=data(*,indexy) ;pokud prepinc1=1, pak histogram ze vsech koeficientu ;jinak pouze z koef. tykajicich se nasi oblasti ;-----------histogram a jeho parametry-----------------------; hist=histogram(korelaceh(0,*),min=-1,max=1,bins=0.01) hodnotyx=findgen(n_elements(hist))*0.01-1 gauss=gaussfit(hodnotyx,hist,koef,nterms=3) mi=koef(1) & sigma=koef(2) horni=mi+3.0*sigma & dolni=mi-3.0*sigma ; ;----------nakresleni korelaci------------------ widget_control,corrdraw2,get_value=okno1 & wset,okno1 indexy1=where((korelace(0,*) lt dolni) or (korelace(0,*) gt horni),count) tvscl,obrazek if count ne 0 then begin korelace1=korelace(*,indexy1) for i=0,count-1 do begin rectangle,korelace1(1,i),korelace1(2,i),krok,krok ;x,y,s,v rectangle,korelace1(3,i),korelace1(4,i),krok,krok endfor endif rectangle1,x,y,krok,krok ;-------------nakresleni histogramu----------- widget_control,histdraw2,get_value=okno2 & wset,okno2 plot,hodnotyx,hist,xrange=[-1,1],xstyle=1,$ title='DISTRIBUTION OF CORRELATION COEFFICIENTS',$ xtitle='correlation coefficient',ytitle='density' oplot,hodnotyx,gauss,color=255 arrow,100,303,140,303,hsize=0,color=255;,/device ;xyouts,100,303,'',/device xyouts,150,300,'gauss fit',/device xyouts,100,280,'mi='+strcompress(string(mi),/remove_all),/device xyouts,100,265,'sigma='+strcompress(string(sigma),/remove_all),/device end 'toggle1':prepinac=1 'toggle2':prepinac=2 'CORECT':begin ;---------smazani predchozi sekce, nahresleni site a textovych poli----------- widget_control,innerbase,/destroy innerbase=widget_base(base,/row) obrazek=readfits(soubory(pocet/2)) s=size(obrazek) crossdraw=widget_draw(innerbase,xsize=s(1),ysize=s(2)) tvscl,obrazek for i=0,xpocet-1 do begin for j=0,ypocet-1 do rectangle,xpocatek+i*krok,ypocatek+j*krok,krok,krok endfor for i=0,xpocet-1 do xyouts,xpocatek+i*krok-5,ypocatek-10,string(i+1),/device,charsize=0.8,alignment=0.5 for i=0,ypocet-1 do xyouts,xpocatek-20,ypocatek+i*krok+5,string(i+1),/device,charsize=0.8,alignment=0.5 corectbase=widget_base(innerbase,/column,xsize=250) corectfield1=cw_field(corectbase,title='coordinates of 1st region (x,y)') corectfield2=cw_field(corectbase,title='coordinates of 2nd region (x,y)') corectfield3=cw_field(corectbase,title='start image (n)') togbase=widget_base(corectbase,/column,/frame,/exclusive) tog1=widget_button(togbase,value='center of brightness',uvalue='toggle1',/no_release) tog2=widget_button(togbase,value='the brightest points',uvalue='toggle2',/no_release) corectbut=widget_button(corectbase,value='CORRELATE',uvalue='CROSSCORECT') widget_control,tog1,/set_button prepinac=1 end 'CROSSCORECT':begin ;zde se provadi korekce podle teziste nebo nejjasnejsiho bodu ;---------parametry oblasti 1 a 2 - x1,y1,x2,y2,krok------------ widget_control,corectfield3,get_value=poradi & poradi=fix(poradi(0)) widget_control,corectfield1,get_value=sourad & sourad=sourad(0) pozice=strpos(sourad,',') & delka=strlen(sourad) xsourad=fix(strmid(sourad,0,pozice)) & ysourad=fix(strmid(sourad,pozice+1,delka-pozice)) x1=(xsourad-1)*krok+xpocatek & y1=(ysourad-1)*krok+ypocatek widget_control,corectfield2,get_value=sourad & sourad=sourad(0) pozice=strpos(sourad,',') & delka=strlen(sourad) xsourad=fix(strmid(sourad,0,pozice)) & ysourad=fix(strmid(sourad,pozice+1,delka-pozice)) x2=(xsourad-1)*krok+xpocatek & y2=(ysourad-1)*krok+ypocatek w=krok & h=krok ;--------vytvoreni poli pro ukladani integralu a teziste--------- rada1=fltarr(pocet);pro ukladani integralu rada2=fltarr(pocet) teziste1=intarr(2,pocet);pro ukladani teziste teziste2=intarr(2,pocet) if prepinac eq 1 then begin ;pokud vybrano teziste iksy=findgen(w+1,h+1) ;pomocne pole pro vypocet teziste iny=findgen(w+1,h+1) for i=0,h do iksy(*,i)=findgen(w+1) for i=0,w do iny(i,*)=findgen(h+1) endif obrazek=readfits(soubory(pocet/2)) s=size(obrazek) animbase=widget_base(/column) animdraw=widget_draw(animbase,xsize=s(1),ysize=s(2)) widget_control,animbase,/realize widget_control,animdraw,get_value=okno & wset,okno ;-------------cyklus - pocitani polohy teziste--------------- for i=0,pocet-1 do begin ;4697 ;0,pocet-1 obr=readfits(soubory(i)) ;nacteme obrazek if i ge poradi then begin ;prepocet souradnic oblasti az od urciteho cisla ;-----poloha teziste nebo nejjasnejsiho bodu pro prvni oblast------ intenzity1=obr(x1:x1+w,y1:y1+h) ;vytrhneme prvni oblast if prepinac eq 1 then begin ;teziste xt1=total(intenzity1*iksy)/total(intenzity1) yt1=total(intenzity1*iny)/total(intenzity1) xt1=round(xt1) & yt1=round(yt1) ;souradnice teziste v oblasti endif else begin ;nejjasnejsi body indexy1=where(intenzity1 eq max(intenzity1)) iksy1=indexy1 mod (w+1) iny1=indexy1/(w+1) xt1=total(iksy1)/n_elements(iksy1) yt1=total(iny1)/n_elements(iny1) endelse teziste1(0,i)=x1+xt1 & teziste1(1,i)=y1+yt1 ;ulozeni polohy teziste x1=x1+xt1-w/2 & y1=y1+yt1-h/2 ;souradnice nove oblasti ;-----poloha teziste nebo nejjasnejsiho bodu pro druhou oblast------ intenzity2=obr(x2:x2+w,y2:y2+h) ;vytrhneme druhou oblast if prepinac eq 1 then begin ;teziste xt2=total(intenzity2*iksy)/total(intenzity2) yt2=total(intenzity2*iny)/total(intenzity2) xt2=round(xt2) & yt2=round(yt2) ;souradnice teziste v oblasti endif else begin ;nejjasnejsi bod indexy2=where(intenzity2 eq max(intenzity2)) iksy2=indexy2 mod (w+1) iny2=indexy2/(w+1) xt2=total(iksy2)/n_elements(iksy2) yt2=total(iny2)/n_elements(iny2) endelse teziste2(0,i)=x2+xt2 & teziste2(1,i)=y2+yt2 ;ulozeni polohy teziste x2=x2+xt2-w/2 & y2=y2+yt2-h/2 ;souradnice nove oblasti endif else begin teziste1(0,i)=x1+w/2 & teziste1(1,i)=y1+h/2 ;pevna poloha teziste a oblasti teziste2(0,i)=x2+w/2 & teziste2(1,i)=y2+h/2 endelse ;---------vypocet integralu v nove oblasti------------ rada1(i)=integruj(obr,x1,y1,w,h)-prumer(i) rada2(i)=integruj(obr,x2,y2,w,h)-prumer(i) ;ulozeni integralu ;---------animace------------- tvscl,obr plots,x1+w/2,y1+h/2,psym=1,color=255,/device plots,x2+w/2,y2+h/2,psym=1,color=255,/device rectangle,x1,y1,w,h rectangle,x2,y2,w,h print,soubory(i) ;vypis endfor ;------------vykresleni grafu---------------- ;----------poloha teziste------------ stary=!p.multi & !p.multi=[0,2,2,0,1] plot,teziste1(0,*),teziste1(1,*),psym=2,title='1st region', $ xtitle='x position',ytitle='y position',$ xrange=[min(teziste1(0,*))-5,max(teziste1(0,*))+5],$ yrange=[min(teziste1(1,*))-5,max(teziste1(1,*))+5],xstyle=1,ystyle=1 plot,teziste2(0,*),teziste2(1,*),psym=2,title='2nd region', $ xtitle='x position',ytitle='y position', $ xrange=[min(teziste2(0,*))-5,max(teziste2(0,*))+5],$ yrange=[min(teziste2(1,*))-5,max(teziste2(1,*))+5],xstyle=1,ystyle=1 !p.multi=[4,2,4,0,1] plot,teziste1(0,*),title='1st region - x position',xrange=[0,pocet],charsize=1.5,$ yrange=[min(teziste1(0,*))-5,max(teziste1(0,*))+5],xstyle=1,ystyle=1 plot,teziste1(1,*),title='1st region - y position',xrange=[0,pocet],charsize=1.5,$ yrange=[min(teziste1(1,*))-5,max(teziste1(1,*))+5],xstyle=1,ystyle=1 plot,teziste2(0,*),title='2nd region - x position',xrange=[0,pocet],charsize=1.5,$ yrange=[min(teziste2(0,*))-5,max(teziste2(0,*))+5],xstyle=1,ystyle=1 plot,teziste2(1,*),title='2nd region - y position',xrange=[0,pocet],charsize=1.5,$ yrange=[min(teziste2(1,*))-5,max(teziste2(1,*))+5],xstyle=1,ystyle=1 !p.multi=stary write_jpeg,'correct2.jpg',tvrd() ;-------------korelace poloh teziste-------------------- tezbase=widget_base(/column) tezdraw=widget_draw(tezbase,xsize=400,ysize=600) widget_control,tezbase,/realize widget_control,tezdraw,get_value=okno & wset,okno lag=findgen(2*pocet-1)-pocet+1 resx=c_correlate(teziste1(0,*),teziste2(0,*),lag) resy=c_correlate(teziste1(1,*),teziste2(1,*),lag) koefx=resx(pocet-1) & koefx=strcompress(string(koefx),/remove_all) koefy=resy(pocet-1) & koefy=strcompress(string(koefy),/remove_all) stary=!p.multi & !p.multi=[0,1,2] plot,lag,resx,title='CORRELATION OF X COORDINATE c0='+koefx,xtitle='lag',ytitle='correlation coefficient',$ yrange=[-1,1],ystyle=1 plot,lag,resy,title='CORRELATION OF y COORDINATE c0='+koefy,xtitle='lag',ytitle='correlation coefficient',$ yrange=[-1,1],ystyle=1 !p.multi=stary write_jpeg,'correct3.jpg',tvrd() ;----------korelace oblasti po korekci korelbase=widget_base(/column) koreldraw=widget_draw(korelbase,xsize=400,ysize=300) widget_control,korelbase,/realize widget_control,koreldraw,get_value=okno & wset,okno lag=findgen(2*pocet-1)-pocet+1 res=c_correlate(rada1,rada2,lag) koef=res(pocet-1) & koef=strcompress(string(koef),/remove_all) plot,lag,res,title='CROSS CORRELATION c0='+koef,xtitle='lag',ytitle='correlation coefficient',$ yrange=[-1,1],ystyle=1 write_jpeg,'correct1.jpg',tvrd() end endcase end pro korel COMMON start,omenu COMMON baze,base,innerbase base=widget_base(xsize=1000,ysize=600,mbar=bar) fmenu=widget_button(bar,value='FILE',/menu) f1=widget_button(fmenu,value='LOAD WORKING DIRECTORY',uvalue='LOAD') f2=widget_button(fmenu,value='QUIT',uvalue='QUIT') omenu=widget_button(bar,value='OPTIONS',/menu) o1=widget_button(omenu,value='whole image',/menu) o11=widget_button(o1,value='count integrals',uvalue='INTEGRALS') o12=widget_button(o1,value='count correlation coefficients',uvalue='KOEFICIENTY') o13=widget_button(o1,value='show regions with high correlation',uvalue='HISTOGRAM') o14=widget_button(o1,value='show regions which correlate with too many regions',uvalue='MNOHO') o2=widget_button(omenu,value='chosen regions',/menu) o21=widget_button(o2,value='auto correlation',uvalue='AUTO') o22=widget_button(o2,value='cross correlation',uvalue='CROSS') o23=widget_button(o2,value='light curve',uvalue='LIGHT') o24=widget_button(o2,value='distribution',uvalue='DIST') o3=widget_button(omenu,value='correction',/menu) o31=widget_button(o3,value='correlation after correction',uvalue='CORECT') innerbase=widget_base(base) widget_control,base,/realize widget_control,omenu,sensitive=0 d=widget_draw(innerbase,xsize=1000,ysize=600) for i=0,999 do arrow,i,0,i,600,hsize=0,color=25100+(i/4.1) tl=3 & vel=15 & ypos=250 xyouts,200,ypos,'K',alignment=0.5,/device,charsize=vel,charthick=tl,color=25340 xyouts,350,ypos,'O',alignment=0.5,/device,charsize=vel,charthick=tl,color=25280 xyouts,500,ypos,'R',alignment=0.5,/device,charsize=vel,charthick=tl,color=25250 xyouts,650,ypos,'E',alignment=0.5,/device,charsize=vel,charthick=tl,color=25180 xyouts,800,ypos,'L',alignment=0.5,/device,charsize=vel,charthick=tl,color=25120 xmanager,'korel',base end