if !d.name eq 'PS' then begin
  device,xsize=18,ysize=16,yoffset=3
  !p.charthick=3 & !p.thick=6 & !x.thick=6 & !y.thick=6
end
;
!p.charsize=2.5
!x.margin=[6.5,1.5]
!y.margin=[3.2,.3]
siz2=1.8
;
;Define the kernel matrix
num=8
kernel_matrix=make_array(num,num,/double,value=0)
t0 = 10.
t_diff = 1.e-0
;
pc_read_ts,obj=ts,/double
t=ts.t
t_index = where (ts.t ge (t0-t_diff) and ts.t le (t0+t_diff))
k100_100 = ts.k100_100
k100_80  = ts.k100_80
k100_60  = ts.k100_60
k100_50  = ts.k100_50
k100_40  = ts.k100_40
k100_30  = ts.k100_30
k100_20  = ts.k100_20
k100_10  = ts.k100_10

k100_100_mean  = mean(k100_100(t_index))
k100_80_mean  = mean(k100_80(t_index))
k100_60_mean  = mean(k100_60(t_index))
k100_50_mean  = mean(k100_50(t_index))
k100_40_mean  = mean(k100_40(t_index))
k100_30_mean  = mean(k100_30(t_index))
k100_20_mean  = mean(k100_20(t_index))
k100_10_mean  = mean(k100_10(t_index))
;
k80_80 = ts.k80_80
k80_60 = ts.k80_60
k80_50 = ts.k80_50
k80_40 = ts.k80_40
k80_30 = ts.k80_30
k80_20 = ts.k80_20
k80_10 = ts.k80_10
k80_80_mean  = mean(k80_80(t_index))
k80_60_mean  = mean(k80_60(t_index))
k80_50_mean  = mean(k80_50(t_index))
k80_40_mean  = mean(k80_40(t_index))
k80_30_mean  = mean(k80_30(t_index))
k80_20_mean  = mean(k80_20(t_index))
k80_10_mean  = mean(k80_10(t_index))
;
k60_60 = ts.k60_60
k60_50 = ts.k60_50
k60_40 = ts.k60_40
k60_30 = ts.k60_30
k60_20 = ts.k60_20
k60_10 = ts.k60_10
k60_60_mean  = mean(k60_60(t_index))
k60_50_mean  = mean(k60_50(t_index))
k60_40_mean  = mean(k60_40(t_index))
k60_30_mean  = mean(k60_30(t_index))
k60_20_mean  = mean(k60_20(t_index))
k60_10_mean  = mean(k60_10(t_index))
k50_50 = ts.k50_50
k50_40 = ts.k50_40
k50_30 = ts.k50_30
k50_20 = ts.k50_20
k50_10 = ts.k50_10
k50_50_mean  = mean(k50_50(t_index))
k50_40_mean  = mean(k50_40(t_index))
k50_30_mean  = mean(k50_30(t_index))
k50_20_mean  = mean(k50_20(t_index))
k50_10_mean  = mean(k50_10(t_index))
k40_40 = ts.k40_40
k40_30 = ts.k40_30
k40_20 = ts.k40_20
k40_10 = ts.k40_10
k40_40_mean  = mean(k40_40(t_index))
k40_30_mean  = mean(k40_30(t_index))
k40_20_mean  = mean(k40_20(t_index))
k40_10_mean  = mean(k40_10(t_index))
;
k30_30 = ts.k30_30
k30_20 = ts.k30_20
k30_10 = ts.k30_10
k30_30_mean  = mean(k30_30(t_index))
k30_20_mean  = mean(k30_20(t_index))
k30_10_mean  = mean(k30_10(t_index))
k20_20 = ts.k20_20
k20_10 = ts.k20_10
k20_20_mean  = mean(k20_20(t_index))
k20_10_mean  = mean(k20_10(t_index))
;
k10_10 = ts.k10_10
k10_10_mean  = mean(k10_10(t_index))
;                         
kernel_matrix[7,0] = k100_100_mean
kernel_matrix[7,1] = k100_80_mean 
kernel_matrix[7,2] = k100_60_mean 
kernel_matrix[7,3] = k100_50_mean 
kernel_matrix[7,4] = k100_40_mean 
kernel_matrix[7,5] = k100_30_mean 
kernel_matrix[7,6] = k100_20_mean 
kernel_matrix[7,7] = k100_10_mean 
;
kernel_matrix[6,1] = k80_80_mean
kernel_matrix[6,2] = k80_60_mean 
kernel_matrix[6,3] = k80_50_mean 
kernel_matrix[6,4] = k80_40_mean 
kernel_matrix[6,5] = k80_30_mean 
kernel_matrix[6,6] = k80_20_mean 
kernel_matrix[6,7] = k80_10_mean 
;
kernel_matrix[5,2] = k60_60_mean
kernel_matrix[5,3] = k60_50_mean 
kernel_matrix[5,4] = k60_40_mean 
kernel_matrix[5,5] = k60_30_mean 
kernel_matrix[5,6] = k60_20_mean 
kernel_matrix[5,7] = k60_10_mean 
;
kernel_matrix[4,3] = k50_50_mean
kernel_matrix[4,4] = k50_40_mean 
kernel_matrix[4,5] = k50_30_mean 
kernel_matrix[4,6] = k50_20_mean 
kernel_matrix[4,7] = k50_10_mean 
;
kernel_matrix[3,4] = k40_40_mean 
kernel_matrix[3,5] = k40_30_mean 
kernel_matrix[3,6] = k40_20_mean 
kernel_matrix[3,7] = k40_10_mean 
;
kernel_matrix[2,5] = k30_30_mean 
kernel_matrix[2,6] = k30_20_mean 
kernel_matrix[2,7] = k30_10_mean 
;
kernel_matrix[1,6] = k20_20_mean 
kernel_matrix[1,7] = k20_10_mean 
;
kernel_matrix[0,7] = k10_10_mean 
;
save,file='kernel_matrix.sav',kernel_matrix
;fo="(8f12.10)"
;print,kernel_matrix,fo=fo
print,kernel_matrix
;radius=[100,80,60,50,40,30,20,10]
radius=[10,20,30,40,50,60,80,100]
minnd=min(kernel_matrix)
maxnd=max(kernel_matrix)
lev=grange(minnd,maxnd,20)
contour,clip(rotate(kernel_matrix,1),minmax(lev)),radius,radius,nlev=30,/fill,ytit='!8r [!7l!6m]',xtit='!8r [!7l!6m]'
;contour,clip(kernel_matrix,minmax(lev)),nlev=30,/fill,ytit='!8r [!7l!6m]',xtit='!8r [!7l!6m]'
;contour,clip(kernel_matrix,minmax(lev)),radius,radius,nlev=30,/fill,ytit='!8r [!7l!6m]',xtit='!8r [!7l!6m]'
COLORBAR,/VERTICAL,NCOLORS=100, POSITION=[.91, 0.5, .94, 0.90],CHARSIZE=1.5,xticks=5,ytickes=5
;
cwd,run
default,texsubdir,''
print,'mv idl.ps ~/tex/nils/particles_coag/fig/'+texsubdir+'kernel_swarm_'+run+'.ps'
;
END