Program Code
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Our data file has 18 coloumns(no_cols) and 4232 rows (no_rows) and ; columns 1(col_time),3(col_field) and 8(col_diff) represents time, intensity ; field and diff respectively. ; Modify variables called no_cols, no_rows, col_time, col_field and col_diff ; according to your data file. Remember in IDL Column/Row number starts ; from 0 (not from 1). ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; This procedure reads the sgt data file (SGT_DATA.DAT) and plots log field ; with respect to diff and time. PRO Start,rawdata=raw_data,no_rows=no_rows,no_cols,col_time,col_field,col_diff openr, LogicalUnit1, 'SGT_DATA.DAT', /get_lun no_cols=18 no_rows=4232 ; give the column numbers of time, field and diff col_time = 1 col_field = 3 col_diff = 8 raw_data = make_array(no_cols,no_rows,/double) readf, LogicalUnit1,raw_data free_lun, LogicalUnit1 ; plot log_field vs diff and time from the original data (Fig. 1) Set_Plot,'x' !p.multi=[0,1,2] plot,raw_data[col_diff,*],alog10(raw_data[col_field,*]),psym=1, $ xtitle='diff', ytitle='log_field',charsize=2.0 plot,raw_data[col_time,*],alog10(raw_data[col_field,*]),psym=1, $ xtitle='time', ytitle='log_field',charsize=2.0 !p.multi = 0 ; now save the plot in the post script format Set_Plot,'ps' !p.multi=[0,1,2] plot,raw_data[col_diff,*],alog10(raw_data[col_field,*]),psym=1, $ xtitle='diff', ytitle='log_field',charsize=2.0 plot,raw_data[col_time,*],alog10(raw_data[col_field,*]),psym=1, $ xtitle='time', ytitle='log_field',charsize=2.0 !p.multi = 0 end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO Value,pbegin=start_time,pend=end_time,min=mindiff,max=maxdiff ; ask for start and end time to be considered read,prompt='Give start time in hour: ',start_time read,prompt='Give end time in hour: ',end_time ; ask for minimum and maximum diff for correction read,prompt='Give minimum diff for correction: ',mindiff read,prompt='Give maximum diff for correction: ',maxdiff end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; This procedure uses the start and end time to create a new data file ; (NEW.DAT) which is a subset of the original file SGT_DATA.DAT PRO Modified,start_time,end_time,mindiff,maxdiff,raw_data,moddata=mod_data, $ n_rows = nrows,no_rows,no_cols,col_time,col_diff,col_field openr, LogicalUnit1, 'SGT_DATA.DAT', /get_lun openw, LogicalUnit2, 'NEW.DAT', /get_lun word='' nrows =0 for i = 0, no_rows-1 do begin readf,LogicalUnit1,word if((raw_data[col_time,i] gt start_time) and (raw_data[col_time,i] lt end_time )) then begin if((raw_data[col_diff,i] ge mindiff) and (raw_data[col_diff,i] lt maxdiff) ) then begin printf,LogicalUnit2,word nrows = nrows+1 endif endif endfor free_lun, LogicalUnit1 free_lun, LogicalUnit2 ; create a new array to hold data from newly created data file mod_data=make_array(no_cols,nrows,/double) openr,logicalunit1,'NEW.DAT', /get_lun readf, logicalunit1,mod_data free_lun, logicalunit1 ; plot new log_field vs diff and time (Fig. 2) Set_Plot,'x' !p.multi=[0,1,2] plot,mod_data[col_diff,*],alog10(mod_data[col_field,*]),psym=1, $ xtitle='diff', ytitle='log_field',charsize=2.0,yticks=3 plot,mod_data[col_time,*],alog10(mod_data[col_field,*]),psym=1, $ xtitle='time', ytitle='log_field',charsize=2.0,yticks=3 !p.multi = 0 ; now save the plot in the post script format Set_Plot,'ps' !p.multi=[0,1,2] plot,mod_data[col_diff,*],alog10(mod_data[col_field,*]),psym=1, $ xtitle='diff', ytitle='log_field',charsize=2.0,yticks=3 plot,mod_data[col_time,*],alog10(mod_data[col_field,*]),psym=1, $ xtitle='time', ytitle='log_field',charsize=2.0,yticks=3 !p.multi = 0 print, ' Press any key to continue... ' ch=get_kbrd(1) print,'' end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO Distribution,mod_data,col_field,col_time,col_diff,nrows log_field = alog10(mod_data[col_field,*]) field_range = max(log_field) - min(log_field) no_bins = 50 field_bin_width = field_range/no_bins local_min = min(log_field) local_max = local_min + field_bin_width field_bin_array = dblarr(no_bins+1) field_bin_array = findgen(no_bins+1)*field_bin_width + min(log_field) field_pts = intarr(no_bins) ii=0.0 for k = 0, no_bins-1 do begin j = 0 for i = 0L, nrows -1 do begin if((log_field(i) ge local_min) and $ (log_field(i) lt local_max)) then j = j+1 if(k eq no_bins-1) then begin if(log_field(i) eq max(log_field))then ii=ii+1 endif endfor field_pts[k] = j+ii local_min = local_max local_max = local_min + field_bin_width endfor ;find those bins which have no of data points greater than 0 newfield_bin_array = where(field_pts ge 1, count) ; create new field_pts with this new_bin_array new_field_pts = field_pts(newfield_bin_array) ; The "field_px" array holds the observed probability distribution P(log_field) field_px = dblarr(count) field_px = new_field_pts/(total(new_field_pts) * field_bin_width) ; create new bin centre newfield_bin_centre = dblarr(count) for i = 0,count -1 do begin k = newfield_bin_array(i) newfield_bin_centre[i] = (field_bin_array[k] + field_bin_array[k+1])/2 endfor ; Plot log(P(log_field)) vs log_field (Fig. 3) Set_Plot,'x' plot_io, newfield_bin_centre, field_px, xtitle='log_field', ytitle='P(log_field)', $ yrange=[.001, 10] ; now save the plot in the post script format Set_Plot,'ps' plot_io, newfield_bin_centre, field_px, xtitle='log_field', ytitle='P(log_field)', $ yrange=[.001, 10] end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO Smoothing, mod_data,nrows,norm_log_field,nvalues, $ col_time,col_field,col_diff,ans,ans_pr tr_diff = transpose(mod_data[col_diff,*]) tr_log_field = transpose(alog10(mod_data[col_field,*])) tr_time = transpose(mod_data[col_time,*]) ; ask for parameter on which smoothing should be done ans_pr='' ans='' print,'Smoothing should be done on diff or time' print,'Answer d for diff or t for time' read,prompt='d or t? :',ans if(ans eq 'd') then begin result = sort(tr_diff) tr_diff = tr_diff(result) tr_log_field = tr_log_field(result) param = tr_diff ans_pr = 'diff' end if (ans eq 't') then begin param = tr_time ans_pr = 'time' end if (( ans ne 'd') and (ans ne 't')) then begin print,'answer should be d or t' stop endif ; create data arrays mu and sigma to hold average and standard deviation of ; log_field mu = dblarr(nrows) sigma = dblarr(nrows) av_param = dblarr(nrows) order = max([10,0.05*nrows]) ; ask for the smoothing (or averaging) interval and create expanded arrays ; for the smoothed (averaged) quantities mu and sigma. ; Note param is an array for either time or diff. read,prompt='Specify the number of data points over which smoothing will be done: ', nvalues if nvalues mod 2 eq 0 then sm_width = nvalues+1 $ else sm_width = nvalues exp_log_field = [ts_fcast(tr_log_field,order,sm_width/2,/backcast), $ tr_log_field, $ ts_fcast(tr_log_field,order,sm_width/2)] exp_param = [ts_fcast(param,order,sm_width/2,/backcast), $ param, $ ts_fcast(param,order,sm_width/2)] for k = 0,nrows -1 do begin mu[k] = total(exp_log_field[k:sm_width+(k-1)])/sm_width av_param[k] = total(exp_param[k:sm_width+(k-1)])/sm_width sigma[k] = stddev(exp_log_field[k:sm_width+(k-1)]) endfor ;plot smoothed data (Fig. 4) Set_Plot,'x' if(ans eq 't') then begin plot,av_param,mu,psym=1,xtitle = 'smoothed time', $ ytitle='smoothed log_field',title='Smoothed Data', $ charsize=2.0 endif else begin plot,av_param,mu,psym=1,xtitle = 'smoothed diff', $ ytitle='smoothed log_field',title='Smoothed Data', $ charsize=2.0 endelse ; now save the plot in the post script format Set_Plot,'ps' if(ans eq 't') then begin plot,av_param,mu,psym=1,xtitle = 'smoothed time', $ ytitle='smoothed log_field',title='Smoothed Data', $ charsize=2.0 endif else begin plot,av_param,mu,psym=1,xtitle = 'smoothed diff', $ ytitle='smoothed log_field',title='Smoothed Data', $ charsize=2.0 endelse print, ' Press any key to continue... ' ch=get_kbrd(1) print,'' ; create normalized data for comparison with SGT norm_log_field = dblarr(nrows) for i = 0, nrows -1 do $ norm_log_field[i] = (tr_log_field[i] - mu[i])/sigma[i] ; plot normalized log_field data vs smoothed time or diff parameter (Fig. 5) Set_Plot,'x' if(ans eq 't') then begin plot,av_param,norm_log_field,psym=1,xtitle='smoothed time ', $ ytitle = 'normalized log_field',title='Normalised Data', $ charsize=2.0 endif else begin plot,av_param,norm_log_field,psym=1,xtitle='smoothed diff', $ ytitle = 'normalized log_field',title='Normalised Data', $ charsize=2.0 endelse ; now save the plot in the post script format Set_Plot,'ps' if(ans eq 't') then begin plot,av_param,norm_log_field,psym=1,xtitle='smoothed time ', $ ytitle = 'normalized log_field',title='Normalised Data', $ charsize=2.0 endif else begin plot,av_param,norm_log_field,psym=1,xtitle='smoothed diff', $ ytitle = 'normalized log_field',title='Normalised Data', $ charsize=2.0 endelse print, ' Press any key to continue... ' ch=get_kbrd(1) print,'' end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO Binwork,norm_log_field,nrows,binwidth=bin_width,numpts=num_pts, $ ptheory=p_theory,new_num_pts,new_bin_array,new_bin_centre, $ new_total_pts,new_no_bins,bin_count ; Now bin the normalized data and compare with SGT using the Chi-Squared ; and Kolmogorov-Smirnov statistical tests print,"X is the normalized log_field variable" read,prompt='bin size in X? ',bin_width bin_max = 5.0 bin_min = -5.0 nx_min = bin_min nx_max = nx_min + bin_width nx_bin = fix((bin_max - bin_min)/bin_width) bin_array = findgen(nx_bin+1)*bin_width - (nx_bin*bin_width/2.0) num_pts = intarr(nx_bin) index = dblarr(nx_bin) ; create bin centre array bin_centre = dblarr(nx_bin) for i = 0,nx_bin -1 do $ bin_centre[i] = (bin_array[i] + bin_array[i+1])/2.0 for k = 0, nx_bin-1 do begin j = 0 for i = 0, nrows -1 do begin if((norm_log_field(i) ge nx_min) and $ (norm_log_field(i) lt nx_max)) then j = j+1 endfor num_pts[k] = j nx_min = nx_max nx_max = nx_max + bin_width endfor print,'...Make a selection about the minimum no of data points in each bin' print,'' print,'...Bins having less than this number will not be considered in the statistical significance test' print,'' print,'...A good choice will be a number greater or equal to 10' print,'' read,prompt='Minimum number of data points in each bin: ',bin_count ; find minimum and maximum bin array depending on bin_count for i=0, nx_bin -1 do begin if (num_pts[i] ge bin_count) then begin min_no = i goto, jumpout1 endif endfor jumpout1: print,'min_no=',min_no for i= nx_bin -1,0,-1 do begin if (num_pts[i] gt 10) then begin max_no = i goto, jumpout2 endif endfor jumpout2: print,'max_no=',max_no ; count the total no of field pts within this restricted domain new_total_pts = 0.0 new_no_bins = 0.0 for i=min_no,max_no do begin new_total_pts = new_total_pts + num_pts[i] new_no_bins = new_no_bins + 1 endfor ; create new bin centre array new_bin_centre = bin_centre[min_no:max_no] new_num_pts = num_pts[min_no:max_no] new_bin_array = bin_array[min_no:(max_no+1)] ; The "index" array holds the observed probability distribution P(x) index = new_num_pts/(total(new_num_pts) * bin_width) ; plot distribution of normalized data Set_Plot,'x' plot,new_bin_centre,index,title='Comparison of observed and predicted P(x)', $ position = [.15,.5,.85,.90],psym=1 xyouts,.65,.80,/Normal,'observed = ++' xyouts,.65,.85,/Normal,'predicted = --' ; plot distribution graph of theory. (Fig. 6) SGT predicts that P(x) is Gaussian ; with zero mean and unit standard deviation p_theory = exp(-0.5*new_bin_centre^2)/2.5066 oplot,new_bin_centre,p_theory ; now save the plot in the post script format Set_Plot,'ps' plot,new_bin_centre,index,title='Comparison of observed and predicted P(x)', $ position = [.15,.5,.85,.90],psym=1 xyouts,.65,.80,/Normal,'observed = ++' xyouts,.65,.85,/Normal,'predicted = --' p_theory = exp(-0.5*new_bin_centre^2)/2.5066 oplot,new_bin_centre,p_theory end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO Chisqtest,p_theory,new_num_pts,bin_width,new_total_pts ; Chi-Squared test from Numerical Recipes by Press et al. expected_val = p_theory*bin_width*new_total_pts observed_val = new_num_pts test_res = xsq_test(observed_val,expected_val,excell=ex_cell, $ obcell = ob_cell) ;draw a box to put all parameters Fig. 6 Set_Plot,'x' box_x_coords = [.15, .15, .85, .85, .15] box_y_coords = [ .2, .45, .45, .2, .2] PlotS, box_x_coords, box_y_coords, /Normal xyouts, .20,.32, /Normal, string(format='("degrees of freedom =",f6.2)', n_elements(ex_cell) -1) xyouts, .20,.30, /Normal, string(format='("chi-square test statistics(z) : ", f12.7)' , test_res[0]) xyouts, .20,.28, /Normal,string(format='("probability of obtaining a value of z or larger: ",f12.7)', test_res[1]) ; now save the plot in the post script format Set_Plot,'ps' box_x_coords = [.15, .15, .85, .85, .15] box_y_coords = [ .2, .45, .45, .2, .2] PlotS, box_x_coords, box_y_coords, /Normal xyouts, .20,.32, /Normal, string(format='("degrees of freedom =",f6.2)', n_ elements(ex_cell) -1) xyouts, .20,.30, /Normal, string(format='("chi-square test statistics(z) : ", f12.7)' , test_res[0]) xyouts, .20,.28, /Normal,string(format='("probability of obtaining a value of z or larger: ",f12.7)', test_res[1]) end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Function Gammln,xx cof = dblarr(6) cof[0] = 76.18009173 cof[1] = -86.50532033 cof[2] = 24.01409822 cof[3] = -1.231739516 cof[4] = .120858003e-2 cof[5] = -.536382e-5 stp = double(2.50662827465) half = double(.5) one = double(1.) fpf = double(5.5) x = double(xx - one) tmp = double(x + fpf) tmp = (x+half)*alog(tmp) -tmp ser = double(one) for i = 0,5 do begin x = x + one ser = ser + cof[i]/x endfor gammln = tmp + alog(stp*ser) return, gammln end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Function Gcf,a,x itmax = 100 eps = 3.e-7 gln = Gammln(a) gold = 0. a0 = 1. a1 = x b0 = 0. b1 = 1. fac = 1. for i = 1,itmax do begin ai = float(i) aia = ai - a a0 = (a1 + a0*aia)*fac b0 = (b1 + b0*aia)*fac aif = ai*fac a1 = x*a0+aif*a1 b1 = x*b0+aif*b1 if(a1 ne 0.) then begin fac = 1./a1 g = b1*fac if( abs((g - gold)/g) lt eps) then goto, jumpout gold = g endif endfor jumpout: gammcf = exp(-x + a * alog(x) -gln)*g return,gammcf end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Function Gser,a,x itmax = 100.0 eps = 3.e-7 gln = Gammln(a) if( x le 0.0) then begin gamser=0. return, gamser endif ap = a sum =1.0/a del = sum for i = 1,itmax do begin ap = ap +1. del = del*x/ap sum = sum + del if(abs(del) lt abs(sum)*eps) then goto, jumpout endfor print,'Atoo large, itmax too small' jumpout:gamser = sum*exp(-x+a*alog(x) - gln) return,gamser end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Function Gammp,a,x if(x lt a+1.) then begin gammp = Gser(a,x) endif else begin gammcf = Gcf(a,x) gammp = 1.- gammcf endelse return, gammp end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Function Gammq,a,x if(x lt a+1.) then begin gamser = Gser(a,x) gammq = 1. - gamser endif else begin gammq = Gcf(a,x) endelse return, gammq end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Function Erf,x if(x lt 0.) then erf = - Gammp( .5, x^2) else $ erf = Gammp( .5, x^2) return,erf end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Function Erfc,x if(x lt 0.) then erfc = 1.0+ Gammp(.5, x^2) else $ erfc = Gammq(.5, x^2) return,erfc end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Function Funcsgt,xx if(xx lt 0.0) then temp = 0.5* Erfc(-0.7071*xx) $ else temp = 0.5 + 0.5* Erf(0.7071*xx) func_sgt = temp return,func_sgt end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; Function Probks,alam eps1 = 0.001 eps2 = 1.e-8 a2 = -2.*alam^2 fac = 2. probks = 0. termbf = 0. for i = 1,100 do begin term = fac*exp(a2*float(i^2)) probks = probks + term if((abs(term) le eps1*termbf) or (abs(term) le eps2*probks)) then begin return, probks endif else begin fac = -fac termbf = abs(term) endelse endfor probks = 1. return,probks end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PRO Ksone,norm_log_etrue_h,nrows,d_sgt=d,prob_sgt=prob,sortdata= sort_norm_data ; Kolmogorov-Smirnov test from Numerical Recipes by Press et al. sort_norm_data = norm_log_etrue_h(sort(norm_log_etrue_h)) en = float(nrows) d = 0. f0 = 0. for i = 1,nrows do begin fn = float(i)/en ff = Funcsgt(sort_norm_data[i-1]) if(abs(f0-ff) ge abs(fn-ff)) then dt = abs(f0-ff) else dt = abs(fn-ff) if(dt gt d) then d=dt f0 = fn endfor prob = Probks(sqrt(nrows)*d) end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; thisDevice = !D.Name Set_Plot,'ps' Device, File='sgt.ps',xsize=18.5,ysize=25.5,xoffset=1.25,yoffset=2.0 Start,rawdata=raw_data,no_rows=no_rows,no_cols,col_time,col_field,col_diff Value,pbegin=start_time,pend=end_time,min=mindiff,max=maxdiff Modified,start_time,end_time,mindiff,maxdiff,raw_data,moddata=mod_data, $ n_rows=nrows,no_rows,no_cols,col_time,col_diff,col_field Distribution,mod_data,col_field,col_time,col_diff,nrows Smoothing,mod_data,nrows,norm_log_field,nvalues, $ col_time,col_field,col_diff,ans,ans_pr Binwork,norm_log_field,nrows,binwidth=bin_width,numpts=num_pts, $ ptheory=p_theory,new_num_pts,new_bin_array,new_bin_centre, $ new_total_pts,new_no_bins,bin_count ; perform Chi Square Test Chisqtest,p_theory,new_num_pts,bin_width,new_total_pts ; perform Kolmogorov-Smirnov Test Ksone,norm_log_field,nrows,d_sgt=d,prob_sgt=prob,sortdata= sort_norm_data Set_Plot,'x' xyouts,.20,.26, /Normal,string(format='("K-S statistic d for sgt = ",f12.7)',d) xyouts,.20,.24, /Normal, string('Probability (d > observed) =',prob) xyouts,.20,.40, /Normal, string(format='("Smoothed over ",A)',ans_pr) xyouts,.20,.38, /Normal,string(format='("Smoothed over",f9.2," data points")',fix(nvalues)) xyouts,.20,.36, /Normal,string('Binsize in X: ',bin_width) xyouts,.20,.34, /Normal,string('No of data sample in each bin: ',bin_count) ;now save the plot in the post script format Set_Plot,'ps' xyouts,.20,.26, /Normal,string(format='("K-S statistic d for sgt = ",f12.7)',d) xyouts,.20,.24, /Normal, string('Probability (d > observed) =',prob) xyouts,.20,.40, /Normal, string(format='("Smoothed over ",A)',ans_pr) xyouts,.20,.38, /Normal,string(format='("Smoothed over",f9.2," data points")',fix(nvalues)) xyouts,.20,.36, /Normal,string('Binsize in X: ',bin_width) xyouts,.20,.34, /Normal,string('No of sample data in each bin: ',bin_count) Device, /Close_File Set_Plot, thisDevice end
