clc
clear all
close all
Desired_Total_Number_Axons = 20 ;
Number_Motor_Axons = 0 * Desired_Total_Number_Axons ;
Number_Sensory_Axons = Desired_Total_Number_Axons - Number_Motor_Axons ;
Number_C = 0 ;
Number_Adelta = 0 ;
Number_Abeta = 20 ;
Total_Number_Axons = Number_Motor_Axons + Number_C + Number_Adelta + Number_Abeta ;
Motor_Axons_radius = 2.85:8 ;
mean_Motor_Axons_radius = mean(Motor_Axons_radius) ;
Sensory_C_radius = 0.1:0.1:0.75 ;
mean_Sensory_C_radius = mean(Sensory_C_radius) ;
Sensory_Adelta_radius = 0.5:0.1:3 ;
mean_Sensory_Adelta_radius = mean(Sensory_Adelta_radius) ;
Sensory_Abeta_radius = 3:0.1:6 ;
mean_Sensory_Abeta_radius = mean(Sensory_Abeta_radius) ;
Sensory_Aalpha_radius = 6:0.1:10 ;
mean_Sensory_Aalpha_radius = mean(Sensory_Aalpha_radius) ;
Motor_Axons_dist = normrnd(mean_Motor_Axons_radius, 0.5 ,[1,Number_Motor_Axons]) ;
for i=1:length(Motor_Axons_dist)
Motor_Axons_dist(1,i) = round(Motor_Axons_dist(1,i),1) ;
end
Motor_Axons_dist ;
Sensory_C_dist = mean_Sensory_C_radius*ones(1,Number_C) ;
Sensory_C_dist ;
figure
plot(1:Number_C, Sensory_C_dist , 'ob')
ylim([0 15])
title('C fibers')
ylabel('radius')
xlabel('#axon')
Sensory_Adelta_dist = abs (normrnd(mean_Sensory_Adelta_radius , 0.5 ,[1,Number_Adelta]) ) ;
for i=1:length(Sensory_Adelta_dist)
Sensory_Adelta_dist(1,i) = round(Sensory_Adelta_dist(1,i),1) ;
end
Sensory_Adelta_dist ;
figure
plot(1:Number_Adelta, Sensory_Adelta_dist , 'ob')
ylim([0 15])
title('Adelta fibers')
ylabel('radius')
Sensory_Abeta_dist = abs (normrnd(mean_Sensory_Abeta_radius , 0.5 ,[1,Number_Abeta]) ) ;
for i=1:length(Sensory_Abeta_dist)
Sensory_Abeta_dist(1,i) = 6 ;
end
Sensory_Abeta_dist ;
figure
plot(1:Number_Abeta, Sensory_Abeta_dist , 'ob')
ylim([0 15])
title('Abeta fibers')
ylabel('radius')
All_Axons(1:Number_Motor_Axons,1) = Motor_Axons_dist ;
All_Axons(Number_Motor_Axons+1:Number_Motor_Axons+Number_C,1) = Sensory_C_dist ;
All_Axons(Number_Motor_Axons+Number_C+1:Number_Motor_Axons+Number_C+Number_Adelta,1) = Sensory_Adelta_dist ;
All_Axons(Number_Motor_Axons+Number_C+Number_Adelta+1:Number_Motor_Axons+Number_C+Number_Adelta+Number_Abeta,1) = Sensory_Abeta_dist ;
size(All_Axons);
GG(1:Number_Motor_Axons,1)=1 ;
GG(Number_Motor_Axons+1:Number_Motor_Axons+Number_C,1)=2 ;
GG(Number_Motor_Axons+Number_C+1:Number_Motor_Axons+Number_C+Number_Adelta,1)=3 ;
GG(Number_Motor_Axons+Number_C+Number_Adelta+1:Number_Motor_Axons+Number_C+Number_Adelta+Number_Abeta,1)=4 ;
shuffled_All_Axons = All_Axons(randperm(length(All_Axons)))' ;
for i=1:length(shuffled_All_Axons)
f = find(All_Axons == shuffled_All_Axons(1,i)) ;
G(i,1)=GG(f(1,1),1) ;
end
G ;
[C,R] = pack_v(shuffled_All_Axons) ;
figure
for i=1:length(R)
if G(i,1)==1
viscircles(C(i,:),R(i,1), 'color', 'b','Linewidth',0.2)
elseif G(i,1)==2
viscircles(C(i,:),R(i,1), 'color', 'r','Linewidth',0.2)
elseif G(i,1)==3
viscircles(C(i,:),R(i,1), 'color', 'g','Linewidth',0.2)
elseif G(i,1)==4
viscircles(C(i,:),R(i,1), 'color', 'k','Linewidth',0.2)
end
hold on
end
axis square
axis square
a = [0:length(C)-1]'; b = num2str(a); c = cellstr(b);
dx = 0.2; dy = 0.1;
text(C(:,1)-dx, C(:,2), c , 'FontSize',10);
for i=2:length(R)
R(i,1) = R(i-1,1) + 0.0000001 ;
end
unique_radius = unique(R) ;
size(unique_radius) ;
for i = 1:length(unique_radius)
x= unique_radius(i,1)*2 ;
p1 = 0.01876 ;
p2 = 0.4787 ;
p3 = 0.1204 ;
axonD = p1*x^2 + p2*x + p3 ;
paraD2 = axonD ;
p1 = 0.006304 ;
p2 = 0.2071 ;
p3 = 0.5339 ;
nodeD = p1*x^2 + p2*x + p3 ;
paraD1 = nodeD ;
p1 = -8.215 ;
p2 = 272.4 ;
p3 = -780.2 ;
deltax = p1*x^2 + p2*x + p3 ;
p1 = -0.0199 ;
p2 = 3.016 ;
p3 = 17.44 ;
paralength2 = p1*x^2 + p2*x + p3 ;
paralength2_segments = 5 ;
p1 = -0.389 ;
p2 = 14.88 ;
p3 = 9.721 ;
nl = round ( p1*x^2 + p2*x + p3 ) ;
paralength1=3 ;
nodelength=1.0 ;
interlength=(deltax-nodelength-(2*paralength1)-(2*paralength2)) ;
internode_segments = 40 ;
parameters(i,1) = axonD ;
parameters(i,2) = nodeD ;
parameters(i,3) = deltax ;
parameters(i,4) = paralength2/paralength2_segments ;
parameters(i,5) = nl ;
parameters(i,6) = interlength/internode_segments ;
end
for ww = 2:20
parameters(ww,1) = parameters(1,1) ;
parameters(ww,2) = parameters(1,2) ;
parameters(ww,3) = parameters(1,3) ;
parameters(ww,4) = parameters(1,4) ;
parameters(ww,5) = parameters(1,5) ;
parameters(ww,6) = parameters(1,6) ;
end
N=zeros(length(C),length(C)) ;
dist=zeros(length(C),length(C)) ;
Distance_edge=zeros(length(C),length(C)) ;
Distance_Centers = pdist2(C, C) ;
for i=1:length(Distance_Centers)
for j=1:length(Distance_Centers)
if i ~= j
Distance_edge(i,j)= Distance_Centers(i,j) - (R(i,1)+R(j,1)) ;
end
end
end
m= 1000 ;
for k=1:length(Distance_Centers)
for i=1:length(Distance_Centers)
if k ~= i
if Distance_Centers(k,i)< m
N(k,i)=1 ;
dist(k,i) = Distance_Centers(k,i) ;
end
end
end
end
for j=1:length(N)
for i=1:length(N)
if N(j,i)==1
N(i,j)=0 ;
dist(i,j) = 0 ;
end
end
end
Total_Length_Axon = 46500 ;
for J=1: length(unique_radius)
z1=0 ;
z2=0 ;
z3=0 ;
z4=0 ;
j=1 ;
i=2 ;
which_radius = J ;
eval(['coordinates_sections_' num2str(J) '(1,1) = 0 ;']);
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+nodelength;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = 0.5;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("node_%d",z1) ;']);
z1=z1+1 ;
i=i+1 ;
j=j+1 ;
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+paralength1;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("MYSA_%d",z2) ;']);
z2=z2+1 ;
i=i+1 ;
j=j+1 ;
for rr=1:paralength2_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+parameters(which_radius,4);']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("FLUT_%d",z3) ;']);
z3=z3+1 ;
i=i+1 ;
j=j+1 ;
end
for r=1:internode_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+(parameters(which_radius,6));']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("STIN_%d",z4) ;']);
z4=z4+1 ;
i=i+1 ;
j=j+1 ;
end
for rr=1:paralength2_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+parameters(which_radius,4);']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("FLUT_%d",z3) ;']);
z3=z3+1 ;
i=i+1 ;
j=j+1 ;
end
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+paralength1;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("MYSA_%d",z2) ;']);
z2=z2+1 ;
i=i+1 ;
j=j+1 ;
Length_one_set(J,1) = eval(['coordinates_sections_' num2str(J) '(1,i-1)']) ;
number_sets(J,1) = round (Total_Length_Axon/Length_one_set(J,1)) ;
for ii=2:number_sets(J,1)
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+nodelength;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("node_%d",z1) ;']);
z1=z1+1 ;
i=i+1 ;
j=j+1 ;
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+paralength1;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("MYSA_%d",z2) ;']);
z2=z2+1 ;
i=i+1 ;
j=j+1 ;
for rr=1:paralength2_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+parameters(which_radius,4);']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("FLUT_%d",z3) ;']);
z3=z3+1 ;
i=i+1 ;
j=j+1 ;
end
for r=1:internode_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+(parameters(which_radius,6));']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("STIN_%d",z4) ;']);
z4=z4+1 ;
i=i+1 ;
j=j+1 ;
end
for rr=1:paralength2_segments
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+parameters(which_radius,4);']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("FLUT_%d",z3) ;']);
z3=z3+1 ;
i=i+1 ;
j=j+1 ;
end
eval(['coordinates_sections_' num2str(J) '(1,i) = coordinates_sections_' num2str(J) '(1,i-1)+paralength1;']);
eval(['mid_coordinates_sections_' num2str(J) '(1,j) = (coordinates_sections_' num2str(J) '(1,i) + coordinates_sections_' num2str(J) '(1,i-1) )/2;']);
eval(['name_sections_' num2str(J) '(1,j) = sprintf("MYSA_%d",z2) ;']);
z2=z2+1 ;
i=i+1 ;
j=j+1 ;
end
Length_All_sets(J,1) = eval(['coordinates_sections_' num2str(J) '(1,i-1)']) ;
Number_of_nodes (J,1) = z1 ;
end
rrr = [67.4678402083890 ;
330.422368924505;
129.465457490639;
219.869784370584;
68.9098714719089;
250.424487623081;
109.396054368700;
272.096904966341;
286.713233306243;
311.231062614663;
187.425304977039;
34.8696932467240;
95.2544189861966;
379.948342384695;
63.3892558911968;
343.539862635652;
223.950453068184;
414.392042116784;
32.5210199613244;
184.154160226586] ;
coordinates_sections_2 = coordinates_sections_2 - rrr(1,1) ;
mid_coordinates_sections_2 = mid_coordinates_sections_2 - rrr(1,1) ;
coordinates_sections_3 = coordinates_sections_3 - rrr(2,1) ;
mid_coordinates_sections_3 = mid_coordinates_sections_3 - rrr(2,1) ;
coordinates_sections_4 = coordinates_sections_4 - rrr(3,1) ;
mid_coordinates_sections_4 = mid_coordinates_sections_4 - rrr(3,1) ;
coordinates_sections_5 = coordinates_sections_5 - rrr(4,1) ;
mid_coordinates_sections_5 = mid_coordinates_sections_5 - rrr(4,1) ;
coordinates_sections_6 = coordinates_sections_6 - rrr(5,1) ;
mid_coordinates_sections_6 = mid_coordinates_sections_6 - rrr(5,1) ;
coordinates_sections_7 = coordinates_sections_7 - rrr(6,1) ;
mid_coordinates_sections_7 = mid_coordinates_sections_7 - rrr(6,1) ;
coordinates_sections_8 = coordinates_sections_8 - rrr(7,1) ;
mid_coordinates_sections_8 = mid_coordinates_sections_8 - rrr(7,1) ;
coordinates_sections_9 = coordinates_sections_9 - rrr(8,1) ;
mid_coordinates_sections_9 = mid_coordinates_sections_9 - rrr(8,1) ;
coordinates_sections_10 = coordinates_sections_10 - rrr(9,1) ;
mid_coordinates_sections_10 = mid_coordinates_sections_10 - rrr(9,1) ;
coordinates_sections_11 = coordinates_sections_11 - rrr(10,1) ;
mid_coordinates_sections_11 = mid_coordinates_sections_11 - rrr(10,1) ;
coordinates_sections_12 = coordinates_sections_12 - rrr(11,1) ;
mid_coordinates_sections_12 = mid_coordinates_sections_12 - rrr(11,1) ;
coordinates_sections_13 = coordinates_sections_13 - rrr(12,1) ;
mid_coordinates_sections_13 = mid_coordinates_sections_13 - rrr(12,1) ;
coordinates_sections_14 = coordinates_sections_14 - rrr(13,1) ;
mid_coordinates_sections_14 = mid_coordinates_sections_14 - rrr(13,1) ;
coordinates_sections_15 = coordinates_sections_15 - rrr(14,1) ;
mid_coordinates_sections_15 = mid_coordinates_sections_15 - rrr(14,1) ;
coordinates_sections_16 = coordinates_sections_16 - rrr(15,1) ;
mid_coordinates_sections_16 = mid_coordinates_sections_16 - rrr(15,1) ;
coordinates_sections_17 = coordinates_sections_17 - rrr(16,1) ;
mid_coordinates_sections_17 = mid_coordinates_sections_17 - rrr(16,1) ;
coordinates_sections_18 = coordinates_sections_18 - rrr(17,1) ;
mid_coordinates_sections_18 = mid_coordinates_sections_18 - rrr(17,1) ;
coordinates_sections_19 = coordinates_sections_19 - rrr(18,1) ;
mid_coordinates_sections_19 = mid_coordinates_sections_19 - rrr(18,1) ;
coordinates_sections_20 = coordinates_sections_20 - rrr(19,1) ;
mid_coordinates_sections_20 = mid_coordinates_sections_20 - rrr(19,1) ;
figure
plot(coordinates_sections_3 ,1 , 'ob')
hold on
plot(coordinates_sections_6 ,1.05 , 'or')
ylim([1 1.5])
compartments_in_a_set = 3 + internode_segments + (2*paralength2_segments) ;
for s=1 :length((unique_radius))
for t=1 :length((unique_radius))
if t> s
clear S_first
clear S_second
clear n nn f1 f2 fz m mm
L1=0 ;
L2=0 ;
p=1 ;
if eval(['coordinates_sections_' num2str(t) '(1,1)']) == eval(['coordinates_sections_' num2str(s) '(1,1)'])
S_first(p,1) = sprintf("node_%d",L1) ;
S_second(p,1) = sprintf("node_%d",L2) ;
p=p+1 ;
L1=L1+1 ;
L2=L2+1 ;
end
for n2=((L2*compartments_in_a_set)+1):compartments_in_a_set:compartments_in_a_set*Number_of_nodes(t,1)
f2 = eval(['mid_coordinates_sections_' num2str(t) '(1,n2)']) ;
for n1=((L1*compartments_in_a_set)+1):compartments_in_a_set:compartments_in_a_set*Number_of_nodes(s,1)
f1 = eval(['mid_coordinates_sections_' num2str(s) '(1,n1)']) ;
if f1 < f2
for m=1:length (eval(['coordinates_sections_' num2str(t)]))-1
if f1 >= eval(['coordinates_sections_' num2str(t) '(1,m)']) && f1 <= eval(['coordinates_sections_' num2str(t) '(1,m+1)'])
S_first(p,1) = sprintf("node_%d",L1) ;
S_second(p,1) = eval(['name_sections_' num2str(t) '(1,m)']) ;
p=p+1 ;
end
end
L1=L1+1 ;
end
end
for mm=1:length (eval(['coordinates_sections_' num2str(s)]))-1
if f2 >= eval(['coordinates_sections_' num2str(s) '(1,mm)']) && f2 <= eval(['coordinates_sections_' num2str(s) '(1,mm+1)'])
S_first(p,1) = eval(['name_sections_' num2str(s) '(1,mm)']) ;
S_second(p,1) = sprintf("node_%d",L2) ;
p=p+1 ;
end
end
L2=L2+1 ;
end
if Number_of_nodes(s,1) > Number_of_nodes(t,1)
for nn=((L1*compartments_in_a_set)+1):compartments_in_a_set: compartments_in_a_set*Number_of_nodes(s,1)
fz=eval(['mid_coordinates_sections_' num2str(s) '(1,nn)']) ;
for m=1:length (eval(['coordinates_sections_' num2str(t)]))-1
if fz >= eval(['coordinates_sections_' num2str(t) '(1,m)']) && fz <= eval(['coordinates_sections_' num2str(t) '(1,m+1)'])
S_first(p,1) = sprintf("node_%d",L1) ;
S_second(p,1) = eval(['name_sections_' num2str(t) '(1,m)']) ;
p=p+1 ;
end
end
L1=L1+1 ;
end
end
eval(['Connect_types_' num2str(s),num2str(t) '(1:length(S_first),1) = S_first ;']);
eval(['Connect_types_' num2str(s),num2str(t) '(length(S_first)+1:length(S_first)+length(S_second),1) = S_second ;']);
clear SAVE
SAVE = char( eval(['Connect_types_' num2str(s),num2str(t)])) ;
filename=['Connect_types_' num2str(s),num2str(t) '.mat'];
save(filename , 'SAVE');
end
end
end
for s=1 :length((unique_radius))
for t=1 :length((unique_radius))
if t== s
clear S_first
clear S_second
clear n nn f1 f2 fz m mm
L1=0 ;
L2=0 ;
p=1 ;
for n=0:Number_of_nodes(s,1)-1
S_first(p,1) = sprintf("node_%d",n) ;
S_second(p,1) = sprintf("node_%d",n) ;
p=p+1 ;
end
eval(['Connect_types_' num2str(s),num2str(t) '(1:length(S_first),1) = S_first ;']);
eval(['Connect_types_' num2str(s),num2str(t) '(length(S_first)+1:length(S_first)+length(S_second),1) = S_second ;']);
clear SAVE
SAVE = char( eval(['Connect_types_' num2str(s),num2str(t)])) ;
filename=['Connect_types_' num2str(s),num2str(t) '.mat'];
save(filename , 'SAVE');
end
end
end
save('R.txt' , 'R' , '-ascii')
save('C.txt' , 'C' , '-ascii')
save('G.txt' , 'G' , '-ascii')
save('parameters.txt' , 'parameters' , '-ascii')
save('Number_of_nodes.txt' , 'Number_of_nodes' , '-ascii')
save('unique_radius.txt' , 'unique_radius' , '-ascii')
save('neighboringAxon.txt' , 'N' , '-ascii')
save('dist.txt' , 'dist' , '-ascii')
save('Distance_edge.txt' , 'Distance_edge' , '-ascii')
z=0 ;
i=2 ;
j=1 ;
number_of_sections = 1000 ;
sectionlength= Total_Length_Axon/number_of_sections ;
for ii=0:number_of_sections-1
Boundary_coordinates_sections(1,1) = 0 ;
Boundary_coordinates_sections(1,i) = Boundary_coordinates_sections(1,i-1) + sectionlength ;
Boundary_mid_coordinates_section(1,j) = (Boundary_coordinates_sections(1,i) + Boundary_coordinates_sections(1,i-1) )/2 ;
Boundary_name_sections(1,j) = sprintf("section_%d",z) ;
i=i+1 ;
z=z+1 ;
j=j+1 ;
end
for t=1 :length((unique_radius))
clear S_first
clear S_second
clear n f1 m
L=0 ;
p=1 ;
for n=1:compartments_in_a_set:compartments_in_a_set*Number_of_nodes(t,1)
f1 = eval(['mid_coordinates_sections_' num2str(t) '(1,n)']) ;
for m=1:length (Boundary_coordinates_sections)-1
if f1 >= Boundary_coordinates_sections(1,m) && f1 <= Boundary_coordinates_sections(1,m+1)
S_first(p,1) = sprintf("node_%d",L) ;
S_second(p,1) = Boundary_name_sections(1,m) ;
p=p+1 ;
end
end
L=L+1 ;
end
eval(['Boundary_to_' num2str(t) '(1:length(S_first),1) = S_first ;']);
eval(['Boundary_to_' num2str(t) '(length(S_first)+1:length(S_first)+length(S_second),1) = S_second ;']);
clear SAVE
SAVE = char( eval(['Boundary_to_' num2str(t)])) ;
filename=['Boundary_to_' num2str(t) '.mat'];
save(filename , 'SAVE');
end
save('Boundary_coordinates.txt' , 'Boundary_coordinates' , '-ascii')
save('Boundary_neighboring.txt' , 'Boundary_neighboring' , '-ascii')
save('Boundary_dist.txt' , 'Boundary_dist' , '-ascii')
for i=1 :length((unique_radius))
for j=1 :length((unique_radius))
if j>i
clear Rg
z=2 ;
L = length(eval(['Connect_types_' num2str(i),num2str(j)]))/2 ;
Rg=zeros(L,1) ;
for Z=2:L-1
CO1 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(Z-1,1)']) );
CO2 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(Z+1,1)']));
g1 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO2)']) ;
Rg(z,1) = (g2-g1)/2 ;
z=z+1 ;
if g2 < g1
msg = 'Error occurred. g2<g1';
error(msg)
end
end
CO1 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(1,1)']) );
CO2 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(2,1)']));
g1 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO2)']) ;
Rg(1,1) = (g2-g1)/2 ;
CO1 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(L-1,1)']) );
CO2 = find( eval(['name_sections_' num2str(i)]) == eval(['Connect_types_' num2str(i),num2str(j) '(L,1)']));
g1 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(i) '(1,CO2)']) ;
Rg(L,1) = (g2-g1)/2 ;
filename=['Rg_' num2str(i),num2str(j) '.txt'];
save(filename , 'Rg' , '-ascii')
end
end
end
for i=1 :length((unique_radius))
for j=2 :length((unique_radius))
if j>i
clear Areas
k=0 ;
L = length(eval(['Connect_types_' num2str(i),num2str(j)]))/2 ;
Areas =zeros(L,1) ;
for Z=1:L
if (eval(['Connect_types_' num2str(i),num2str(j) '(Z,1)']) == sprintf("node_%d",k) )
Areas(Z,1) = pi*parameters(i,2) ;
k=k+1 ;
else
Areas(Z,1) = pi*parameters(j,2) ;
end
end
filename=['Areas_' num2str(i),num2str(j) '.txt'];
save(filename , 'Areas' , '-ascii')
end
end
end
for j=1:length((unique_radius))
clear Rg
z=2 ;
L = length(eval(['Boundary_to_' num2str(j)]))/2 ;
Rg=zeros(L,1) ;
for Z=2:L-1
CO1 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(Z-1,1)']) );
CO2 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(Z+1,1)']) );
g1 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO2)']) ;
Rg(z,1) = (g2-g1)/2 ;
z=z+1 ;
end
CO1 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(1,1)']) );
CO2 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(2,1)']));
g1 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO2)']) ;
Rg(1,1) = (g2-g1)/2 ;
CO1 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(L-1,1)']) );
CO2 = find( eval(['name_sections_' num2str(j)]) == eval(['Boundary_to_' num2str(j) '(L,1)']));
g1 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO1)']) ;
g2 = eval(['mid_coordinates_sections_' num2str(j) '(1,CO2)']) ;
Rg(L,1) = (g2-g1)/2 ;
filename=['Boundary_Rg_' num2str(j) '.txt'];
save(filename , 'Rg' , '-ascii')
end