-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathstability_analysis.m
669 lines (616 loc) · 30.4 KB
/
stability_analysis.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
clear all
close all
%% Information abput the code
%This code can perform the traditional kinematic analysis based onto the
%sterographic projection of orientation/attitude data of rock disconinuity.
%It needs about dip, dip direction and set of disconinuities planes and
%intersection. To calculate the intersection in the appropriate way,
%please use the 'intersection_calculator.m' code.
%This works include some sciprts presented of Markovaara-Koivisto (2012)
%for what concerns the sterographic projection
%% Define slope orientation dipdir/dip and friction and lateral limit values
%values dipidir/dip of the slope 'parete3 Ormea' are 300/69
disp('REMEBER TO DEFINE SLOPE ATTUTIDE AND OTHER PARAMETERS FROM LINE 11 TO 24')
slope_min_dip= 75;%min. dip of the slope
slope_max_dip= 75;%max. dip of the slope
dipstep=1;%if you want to analyze only one step of dip, min and max dip must be equal
slope_min_dipdir=280;%min. dip direction of the slope
slope_max_dipdir=280;%max. dip direction of the slope
dipdirstep=1;%if you want to analyze only one step of dip direction, min and max dip direction must be equal
%Define friction angle
friction = 30;
%Define lateral limits
latlimit = 20;
plot=1; %specify 1 or 0 if you want or not the plots
%% Select data to load
%plane data with set
[filename, pathname] = uigetfile({'*.xlsx', 'Select the XLSX file of discontinuity plane with set interpretation'},'Select XLSX file of discontinuity plane with set interpretation',...
'F:\Menegoni\Ormea\Nuovi voli giugno 2017\Albris_SenseFly\parete\obj\measures');
% pathname='F:\Menegoni\Ormea\Nuovi voli giugno 2017\Albris_SenseFly\parete\obj\measures\Fit_Matlab\interpretationset';
% filename='Group.xlsx';
data=readtable(fullfile(pathname,filename));
% uiwait(msgbox(['Data of discontinuity Intersections have been select' newline 'SELECT INTERSECTION DISCONTINUITY DATA']));
%Intersection data
[Int_filename, Int_pathname] = uigetfile({'*.xlsx', 'Select the XLSX file of the discontinutiy intersections'},'Select XLSX file of the discontinutiy intersections',...
'F:\Menegoni\Ormea\Nuovi voli giugno 2017\Albris_SenseFly\parete\obj\measures');
% Int_filename='Intersection_Fit_Group..xlsx';
% Int_pathname='F:\Menegoni\Ormea\Nuovi voli giugno 2017\Albris_SenseFly\parete\obj\measures\Fit_Matlab\intersections_Group';
Int_data=readtable(fullfile(Int_pathname,Int_filename));
%Selection of folder in which images will be saved
if plot == 1
uiwait(msgbox('Select the path in which images will be saved'));
pathFig=uigetdir([...
'F:\Menegoni\Ormea\Nuovi voli giugno 2017\Albris_SenseFly\parete\obj\measures\Fit_Matlab']);
end
%
%% Read data and write Dip, Dip Direction, Set, Pole Dip Direction and Dip vectors
%Discontinuity planes
num_disc=numel(data)/3;
Dip=data.Dip(:);%read Dip value of discontinuities
DipDir=data.DipDirection(:);%read Dip Direction value of discontinuities
Set=data.Set(:);%read Set value of discontinuities (random values are defined as NaN as default)
Set(isnan(Set))=0;%change Set value of random discontinuities from NaN to 0
Set_name=unique(Set);
nSet=numel(unique(Set));
for i= 1 : num_disc %calculate discontinuity line pole (dip and dip direction)
%this is an trick to plot plunge and trend line vector
pole_Dip(i,1) = 90 - Dip(i);%discontinuity pole dip
if DipDir(i) < 180%discontinuity pole dip direction
pole_DipDir(i,1) = DipDir(i) + 180;
else
pole_DipDir(i,1) = DipDir(i) - 180;
end
end
%Discontinuity Intersection
num_Int=numel(Int_data.Trend);
Trend=Int_data.Trend;
Plunge=Int_data.Plunge;
Int_Sets(:,1)=Int_data.Set_i;
Int_Sets(:,2)=Int_data.Set_j;
for i= 1 : num_Int %calculate discontinuity line pole (dip and dip direction)
%this is an trick to plot plunge and trend line vector
pole_Plunge(i,1) = 90 - Plunge(i);%discontinuity pole dip
if Trend(i) < 180%discontinuity pole dip direction
pole_Trend(i,1) = Trend(i) + 180;
else
pole_Trend(i,1) = Trend(i) - 180;
end
end
%Plot discontinuity plane poles and discontinuity intersections lower hemisphere projections
if plot==1
Stereogram(DipDir, Dip, num_disc, 'Discontinutiy plane poles', 0, Set )
Stereogram(pole_Trend, pole_Plunge, num_Int, 'Discontinuity intersection', 0, Trend )
end
disp('DipDir Dip %PS %WS %FT %DT %OT')
%% Kynematic analysis
for dipdirloop = slope_min_dipdir: dipdirstep : slope_max_dipdir
for diploop = slope_min_dip: dipstep : slope_max_dip
clearvars planeslope polesslope strikeslope sxlatlimit dxlatlimit
planeslope = zeros(1,2);
planeslope = [dipdirloop, diploop];
%Setting properly planeslope and poleslope
if planeslope(1,1) < 180%calculate poles dip and dip direction of the slope
poleslope = [planeslope(1,1)+180, 90 - planeslope(1,2)];
else
poleslope = [planeslope(1,1)-180, 90 - planeslope(1,2)];
end
if planeslope(1,1)<90%calculate strike of the slope
strikeslope = 360 + planeslope(1,1) - 90;
else
strikeslope = planeslope(1,1) - 90;
end
if planeslope(1,1) < latlimit
sxlatlimit = 360 + planeslope(1,1) - latlimit;%sxlatlimit=laterl limit at left of dip direction of slope
dxlatlimit = planeslope(1,1) + latlimit;%dxlatlimit=laterl limit at rigth of dip direction of slope
elseif planeslope(1,1) >= 360 - latlimit
sxlatlimit = planeslope(1,1) - latlimit;
dxlatlimit = planeslope(1,1) + latlimit - 360;
else
sxlatlimit = planeslope(1,1) - latlimit;
dxlatlimit = planeslope(1,1) + latlimit;
end
%% Planare sliding (PS) kinematic analysis
%A plane could act as a sliding plane if is dip-slope and dipping with an
%angle lower than apparent angle of the slope along dipdirection of the plane
%and higher than frction angle.
app_angle=zeros(num_disc,1);
PlanarSliding=zeros(num_disc, 1);
for i = 1 : num_disc % calculate if discontinuity plane could be a sliding plane
if (DipDir(i) > strikeslope && DipDir(i) < strikeslope + 180)%Calculate apparent angle of the slope along Dip Direction of the discontinuity
app_angle(i,1) = atand(tand(planeslope(1,2)) * cosd(DipDir(i)- planeslope(1,1)));
elseif (strikeslope+180> 360 && DipDir(i) < strikeslope - 180)
app_angle(i,1) = atand(tand(planeslope(1,2)) * cosd(360 + DipDir(i)- planeslope(1,1)));
else
app_angle(i,1)=NaN;
end
if Dip(i) > friction && Dip(i) < app_angle(i)
if (DipDir(i)>sxlatlimit && DipDir(i)<dxlatlimit) || (sxlatlimit >dxlatlimit && DipDir(i)<sxlatlimit && DipDir(i)<dxlatlimit)|| (sxlatlimit >dxlatlimit && DipDir(i)>sxlatlimit && DipDir(i)>dxlatlimit)
%1st case if sx is < dx latlimit, 2nd and 3rd when dx < sx
%lat limit with DipDir of discontinuity >340 && <360 (3rdcase) and
%DipDir >0 && <20 (2nd case).
PlanarSliding(i) = 2;%value assigned to critical discontinuity inside the lateral limits
else
PlanarSliding(i) = 1;%value assigned to crtitical discontinuity outside the lateral limits
end
else
PlanarSliding(i) = NaN;
end
end
PS_Hcrit=zeros;
PS_Lcrit=zeros;
%Calculate and write statistic
PS_Hcrit=sum(PlanarSliding(:) == 2);
PS_Lcrit=sum(PlanarSliding(:) == 1);
if plot==1
Stereogram(DipDir, Dip, num_disc, 'PlanarSliding', 0, PlanarSliding )
hold on
title(['PlanarSliding, N = ', num2str(num_disc) newline 'critical with lateral limits (red) =', num2str(PS_Hcrit) newline 'critical without lateral limits (colored) =', num2str(PS_Hcrit + PS_Lcrit)])
hold off
end
saveas(gcf, fullfile(pathFig,'PlanarSliding.fig'))
saveas(gcf, fullfile(pathFig,'PlanarSliding.png'))
%remember to write a code to export the results.
%% Flexural Toppling (FT) kinematic analysis
slopelimit=[planeslope(1,1), planeslope(1,2)-friction];
pole_app_angle=zeros(num_disc,1);
FlexuralToppling = zeros(num_disc,1);
for i = 1 : num_disc
%Calculate apparent angle of the slope along Dip Direction of the discontinuity
if (pole_DipDir(i)>sxlatlimit && pole_DipDir(i)<dxlatlimit) || (sxlatlimit >dxlatlimit && pole_DipDir(i)<sxlatlimit && pole_DipDir(i)<dxlatlimit)|| (sxlatlimit >dxlatlimit && pole_DipDir(i)>sxlatlimit && pole_DipDir(i)>dxlatlimit)
%if cylcle for Pole DipDirections included in the lateral limit
if (pole_DipDir(i) > strikeslope && pole_DipDir(i) < strikeslope + 180)
pole_app_angle(i,1) = atand(tand(slopelimit(1,2)) * cosd(pole_DipDir(i)- slopelimit(1,1)));
elseif (strikeslope+180> 360 && pole_DipDir(i) < strikeslope - 180)
pole_app_angle(i,1) = atand(tand(slopelimit(1,2)) * cosd(360 + pole_DipDir(i)- slopelimit(1,1)));
else
pole_app_angle(i,1)=NaN;
end
end
%Calculate if discontinuity could be affected by flexural toppling
if pole_Dip(i) < pole_app_angle(i)
FlexuralToppling(i) = 2;
else
FlexuralToppling(i) = NaN;
end
end
%Calculate and write statistic
FT_crit=zeros;
FT_crit=sum(FlexuralToppling(:) == 2);
if plot==1
Stereogram(DipDir, Dip, num_disc, 'FlexuralToppling', 0, FlexuralToppling )
hold on
title(['FlexuralToppling, N = ', num2str(num_disc) newline 'critical with lateral limits =', num2str(FT_crit)])
hold off
end
saveas(gcf, fullfile(pathFig,'FlexuralToppling.fig'))
saveas(gcf, fullfile(pathFig,'FlexuralToppling.png'))
%% WedgeSliding (WS) kinematic analysis
Int_app_angle=zeros(num_Int,1);
Int_app_angle_friction=zeros(num_Int,1);
WedgeSliding=zeros(num_Int, 1);
for i = 1 : num_Int % calculate if discontinuity plane could be a sliding plane
if (Trend(i) > strikeslope && Trend(i) < strikeslope + 180)%Calculate apparent angle of the slope along Dip Direction of the discontinuity
Int_app_angle(i,1) = atand(tand(planeslope(1,2)) * cosd(Trend(i)- planeslope(1,1)));
Int_app_angle_friction(i,1) = atand(tand(friction) * cosd(Trend(i)- planeslope(1,1)));
elseif (strikeslope+180> 360 && Trend(i) < strikeslope - 180)
Int_app_angle(i,1) = atand(tand(planeslope(1,2)) * cosd(360 + Trend(i)- planeslope(1,1)));
Int_app_angle_friction(i,1) = atand(tand(friction) * cosd(Trend(i)- planeslope(1,1)));
else
Int_app_angle(i,1)=NaN;
Int_app_angle_friction(i,1)=NaN;
end
if Plunge(i) < Int_app_angle(i)
if Plunge(i)>friction
WedgeSliding(i) = 2;%value assigned to critical discontinuity inside the lateral limits
elseif Plunge(i)<friction && Plunge(i)>Int_app_angle_friction(i)
WedgeSliding(i) = 1;%value assigned to crtitical discontinuity outside the lateral limits
else
WedgeSliding(i) = NaN;
end
else
WedgeSliding(i) = NaN;
end
end
%Calculate and write statistic
WS_Hcrit=zeros;
WS_Lcrit=zeros;
WS_Hcrit=sum(WedgeSliding(:) == 2);
WS_Lcrit=sum(WedgeSliding(:) == 1);
if plot==1
Stereogram(pole_Trend, pole_Plunge, num_Int, 'WedgeSliding', 1, WedgeSliding )
hold on
title(['WedgeSliding, N = ', num2str(num_Int) newline ...
'wedge that could slide on both the parental discontinuies (red) =', num2str(WS_Hcrit) newline ...
'critical without lateral limits (pink) =', num2str(WS_Lcrit)])
hold off
end
saveas(gcf, fullfile(pathFig,'WedgeSliding.fig'))
saveas(gcf, fullfile(pathFig,'WedgeSliding.png'))
%% DirectToppling (DT) and ObliqueToppling (OT) (with BasalPlane, BP) kinematic analysis
DirectToppling=zeros(num_Int,1);
ObliqueToppling=zeros(num_Int,1);
BasalPlane=zeros(num_disc,1);
upperTopplingLimit = 90 - planeslope(1,2);%angle between slope face and intersection
%(ho qualche dubbio sul considerare questo limite costante e non variabile
%in base all'inclinazione apparente della parete lungo la direzione
%d'immersione dell'intersezione.
if sxlatlimit(1,1) < 180%calculate the opposite of the lateral limit (sx)
oppsxlatlimit = sxlatlimit+180;
else
oppsxlatlimit = sxlatlimit-180;
end
if dxlatlimit(1,1) < 180%calculate the opposite of the lateral limits (dx)
oppsxlatlimit = dxlatlimit+180;
else
oppdxlatlimit = dxlatlimit-180;
end
for i = 1 : num_Int %Discontinuity Intersection
if (Trend(i) > oppsxlatlimit && Trend(i) < oppdxlatlimit) || ...
(oppdxlatlimit < oppsxlatlimit && Trend(i) > oppsxlatlimit && Trend(i) > oppdxlatlimit) || ...
(oppdxlatlimit < oppsxlatlimit && Trend(i) < oppsxlatlimit && Trend(i) < oppdxlatlimit)
if Plunge(i) > upperTopplingLimit && Plunge(i) < (90 - friction)
DirectToppling(i)=2;
elseif Plunge(i) > (90 - friction)
DirectToppling(i)=1;
else
DirectToppling(i)=NaN;
end
elseif (Plunge(i) > (90 - friction) && strikeslope <= 180 && planeslope(1,1)<= 180 && (Trend(i) < strikeslope || Trend(i)> (strikeslope+180))) || ...
(Plunge(i) > (90 - friction) && strikeslope>= 180 && (planeslope(1,1)<= 360 || planeslope(1,1)<= 0) && Trend(i)>(strikeslope-180) && Trend(i)<strikeslope)
ObliqueToppling(i)=1;
DirectToppling(i)=NaN;
else
DirectToppling(i)=NaN;
ObliqueToppling(i)=NaN;
end
end
for i = 1 : num_disc %Discontinuity Planes
if (pole_DipDir(i) > oppsxlatlimit && pole_DipDir(i) < oppdxlatlimit) || ...
(oppdxlatlimit < oppsxlatlimit && pole_DipDir(i) > oppsxlatlimit && pole_DipDir(i) > oppdxlatlimit) || ...
(oppdxlatlimit < oppsxlatlimit && pole_DipDir(i) < oppsxlatlimit && pole_DipDir(i) < oppdxlatlimit)
if pole_Dip(i) > upperTopplingLimit && pole_Dip(i) < (90 - friction)
BasalPlane(i) = 2;
elseif pole_Dip(i)>(90-friction)
BasalPlane(i) = 1;
else
BasalPlane(i) = NaN;
end
elseif (pole_Dip(i) > (90 - friction) && strikeslope <= 180 && planeslope(1,1)<= 180 && (pole_DipDir(i) < strikeslope || pole_DipDir(i)> (strikeslope+180))) || ...
(pole_Dip(i) > (90 - friction) && strikeslope>= 180 && (planeslope(1,1)<= 360 || planeslope(1,1)<= 0) && pole_DipDir(i)>(strikeslope-180) && pole_DipDir(i)<strikeslope)
BasalPlane(i) = 3;
else
BasalPlane(i)=NaN;
end
end
%DirectToppling
DT_Hcrit=zeros;
DT_Lcrit=zeros;
DT_Hcrit=sum(DirectToppling(:) == 2);
DT_Lcrit=sum(DirectToppling(:) == 1);
if plot==1
Stereogram(pole_Trend, pole_Plunge, num_Int, 'DirectToppling', 1, DirectToppling )
hold on
title(['DirectToppling, N = ', num2str(num_Int) newline ...
'wedge that could be affected by direct toppling and with a sliding release mode (red) =', num2str(DT_Hcrit) newline ...
'wedge that could be affected by direct toppling and without a sliding release mode (pink) =', num2str(DT_Lcrit)])
hold off
end
saveas(gcf, fullfile(pathFig,'DirectToppling.fig'))
saveas(gcf, fullfile(pathFig,'DirectToppling.png'))
%ObliqueToppling
OT_crit=sum(ObliqueToppling(:) == 1);
if plot==1
Stereogram(pole_Trend, pole_Plunge, num_Int, 'ObliqueToppling', 1, ObliqueToppling )
hold on
title(['ObliqueToppling, N = ', num2str(num_Int) newline ...
'wedge that could be affected by oblique toppling and without a sliding basal plane (pink) =', num2str(OT_crit)])
hold off
end
saveas(gcf, fullfile(pathFig,'ObliqueToppling.fig'))
saveas(gcf, fullfile(pathFig,'ObliqueToppling.png'))
%BasalPlane
BP_Hcrit=zeros;
BP_Mcrit=zeros;
BP_Lcrit=zeros;
BP_Hcrit=sum(DirectToppling(:) == 2);
BP_Mcrit=sum(DirectToppling(:) == 1);
BP_Lcrit=sum(BasalPlane(:) == 3);
if plot==1
Stereogram(DipDir, Dip, num_disc, 'BasalPlane', 0, BasalPlane )
hold on
title(['BasalPlane, N = ', num2str(num_Int) newline ...
'BasalPlane that could be need to activate direct toppling with sliding release mode(red) =', num2str(BP_Hcrit) newline...
'BasalPlane that could be need to activate direct toppling without sliding release mode(yellow) =', num2str(BP_Mcrit) newline...
'BasalPlane that could be need to activate oblique toppling without sliding release mode(pink) =', num2str(BP_Lcrit)])
hold off
end
saveas(gcf, fullfile(pathFig,'BasalPlane.fig'))
saveas(gcf, fullfile(pathFig,'BasalPlane.png'))
%% Calculation statistics and display them
%PlanarSliding-Statistic
Set_PS=zeros(nSet+1,5);
critic_PS_Hc=zeros;
critic_PS_Lc=zeros;
for k = 1 : nSet
Hc=0;
Lc=0;
Set_PS(k,1) = Set_name(k);%set number value
Set_PS(k,2) = sum(Set(:) == (Set_name(k)));%total number of discontinuity of set
for i = 1 : num_disc%value(%) of critical discontinuity
if Set(i)== Set_name(k) && PlanarSliding(i) == 2
Hc = Hc + 1;
critic_PS_Hc=critic_PS_Hc+1;
elseif Set(i)== Set_name(k) == 1 && PlanarSliding(i) == 1
Lc=Lc+1;
critic_PS_Lc=critic_PS_Lc+1;
end
end
Set_PS(k,4) =Hc ;
Set_PS(k,5) =Lc ;
Set_PS(k,3) =Set_PS(k,4) + Set_PS(k,5);
end
for i = 1 : numel(Set_PS(1,:))
if i == 1
Set_PS(nSet+1,i)=NaN;
elseif i==2
Set_PS(nSet+1,i)=sum(Set_PS(:,i));
elseif i==3
Set_PS(nSet+1,i)=(critic_PS_Hc+critic_PS_Lc);
elseif i==4
Set_PS(nSet+1,i)=(critic_PS_Hc);
elseif i==5
Set_PS(nSet+1,i)=(critic_PS_Lc);
end
end
clearvars PS_T
PS_T=table(Set_PS(:,1),Set_PS(:,2),Set_PS(:,3),Set_PS(:,4),Set_PS(:,5));
PS_T.Properties.VariableNames{'Var1'} = 'Set';
PS_T.Properties.VariableNames{'Var2'} = 'TotDiscontinuities';
PS_T.Properties.VariableNames{'Var3'} = 'TotCritical';
PS_T.Properties.VariableNames{'Var4'} = 'HighCritical';
PS_T.Properties.VariableNames{'Var5'} = 'LowCritical';
%WedgeSliding-Statistic
Set_WS=zeros(1,6);
index_Set_WS=zeros;
tot_int_sets=zeros;
critic_WS_Hc=zeros;
critic_WS_Lc=zeros;
for k = 1 : nSet
for l = k : nSet
Set_WS(index_Set_WS+1,1) = Set_name(k);
Set_WS(index_Set_WS+1,2) = Set_name(l);
index_Set_WS = numel(Set_WS(:,1));
tot_int_sets = 0;
Hc = 0;
Lc = 0;
for i = 1 : num_Int
if (Int_Sets(i,1) == Set_name(k) && Int_Sets(i,2) == Set_name(l) )||(Int_Sets(i,2) == Set_name(k) && Int_Sets(i,1) == Set_name(l))
tot_int_sets = tot_int_sets + 1;
if WedgeSliding(i) == 2
Hc = Hc + 1;
critic_WS_Hc=critic_WS_Hc+1;
elseif WedgeSliding(i) == 1
Lc = Lc +1;
critic_WS_Lc=critic_WS_Lc+1;
end
end
end
Set_WS(index_Set_WS,3) = tot_int_sets;
Set_WS(index_Set_WS,4) = (Hc + Lc);
Set_WS(index_Set_WS,5) = Hc;
Set_WS(index_Set_WS,6) = Lc;
end
end
for i = 1 : numel(Set_WS(1,:))
if i ==1 || i==2
Set_WS(index_Set_WS+1,i)=NaN;
elseif i==3
Set_WS(index_Set_WS+1,i)=sum(Set_WS(:,i));
elseif i==4
Set_WS(index_Set_WS+1,i)=(critic_WS_Hc+critic_WS_Lc);
elseif i==5
Set_WS(index_Set_WS+1,i)=(critic_WS_Hc);
elseif i==6
Set_WS(index_Set_WS+1,i)=(critic_WS_Lc);
end
end
clearvars WS_T
WS_T=table(Set_WS(:,1),Set_WS(:,2),Set_WS(:,3),Set_WS(:,4),Set_WS(:,5),Set_WS(:,6));
WS_T.Properties.VariableNames{'Var1'} = 'Set_i';
WS_T.Properties.VariableNames{'Var2'} = 'Set_j';
WS_T.Properties.VariableNames{'Var3'} = 'TotIntersections';
WS_T.Properties.VariableNames{'Var4'} = 'TotCritical';
WS_T.Properties.VariableNames{'Var5'} = 'HighCritical';
WS_T.Properties.VariableNames{'Var6'} = 'LowCritical';
%FlexuraToppling-Statistic
Set_FT=zeros(nSet+1,5);
critic_FT=zeros;
for k = 1 : nSet
Hc=0;
Set_FT(k,1) = Set_name(k);%set number value
Set_FT(k,2) = sum(Set(:) == (Set_name(k)));%total number of discontinuity of set
for i = 1 : num_disc%value(%) of critical discontinuity
if Set(i)== Set_name(k) && FlexuralToppling(i) == 2
Hc = Hc + 1;
critic_FT=critic_FT+1;
end
end
Set_FT(k,3) =Hc ;
end
for i = 1 : numel(Set_FT(1,:))
if i == 1
Set_FT(nSet+1,i)=NaN;
elseif i==2
Set_FT(nSet+1,i)=sum(Set_FT(:,i));
else
Set_FT(nSet+1,i)=critic_FT;
end
end
clearvars FT_T
FT_T=table(Set_FT(:,1),Set_FT(:,2),Set_FT(:,3));
FT_T.Properties.VariableNames{'Var1'} = 'Set';
FT_T.Properties.VariableNames{'Var2'} = 'TotDiscontinuities';
FT_T.Properties.VariableNames{'Var3'} = 'TotCritical';
%DirectToppling-Statistic
Set_DT=zeros(1,6);
index_Set_DT=zeros;
tot_int_sets=zeros;
critic_DT=zeros;
critic_DT_Hc=zeros;
critic_DT_Lc=zeros;
for k = 1 : nSet
for l = k : nSet
Set_DT(index_Set_DT+1,1) = Set_name(k);
Set_DT(index_Set_DT+1,2) = Set_name(l);
index_Set_DT = numel(Set_DT(:,1));
tot_int_sets = 0;
Hc = 0;
Lc = 0;
for i = 1 : num_Int
if (Int_Sets(i,1) == Set_name(k) && Int_Sets(i,2) == Set_name(l) )||(Int_Sets(i,2) == Set_name(k) && Int_Sets(i,1) == Set_name(l))
tot_int_sets = tot_int_sets + 1;
if DirectToppling(i) == 2
Hc = Hc + 1;
critic_DT_Hc=critic_DT_Hc+1;
elseif DirectToppling(i) == 1
Lc = Lc +1;
critic_DT_Lc=critic_DT_Lc+1;
end
end
end
Set_DT(index_Set_DT,3) = tot_int_sets;
Set_DT(index_Set_DT,4) = (Hc + Lc);
Set_DT(index_Set_DT,5) = Hc;
Set_DT(index_Set_DT,6) = Lc;
end
end
for i = 1 : numel(Set_DT(1,:))
if i ==1 || i==2
Set_DT(index_Set_DT+1,i)=NaN;
elseif i==3
Set_DT(index_Set_DT+1,i)=sum(Set_DT(:,i));
elseif i==4
Set_DT(index_Set_DT+1,i)=(critic_DT_Hc+critic_DT_Lc);
elseif i==5
Set_DT(index_Set_DT+1,i)=critic_DT_Hc;
elseif i==6
Set_DT(index_Set_DT+1,i)=critic_DT_Lc;
end
end
clearvars DT_T
DT_T=table(Set_DT(:,1),Set_DT(:,2),Set_DT(:,3),Set_DT(:,4),Set_DT(:,5),Set_DT(:,6));
DT_T.Properties.VariableNames{'Var1'} = 'Set_i';
DT_T.Properties.VariableNames{'Var2'} = 'Set_j';
DT_T.Properties.VariableNames{'Var3'} = 'TotIntersections';
DT_T.Properties.VariableNames{'Var4'} = 'TotCritical';
DT_T.Properties.VariableNames{'Var5'} = 'SlidingFailiureMode';
DT_T.Properties.VariableNames{'Var6'} = 'NoSlidingFailureMode';
%ObliqueToppling-Statistic
Set_OT=zeros(1,4);
index_Set_OT=zeros;
tot_int_sets=zeros;
critic_OT=zeros;
for k = 1 : nSet
for l = k : nSet
Set_OT(index_Set_OT+1,1) = Set_name(k);
Set_OT(index_Set_OT+1,2) = Set_name(l);
index_Set_OT = numel(Set_OT(:,1));
tot_int_sets = 0;
Hc = 0;
for i = 1 : num_Int
if (Int_Sets(i,1) == Set_name(k) && Int_Sets(i,2) == Set_name(l) )||(Int_Sets(i,2) == Set_name(k) && Int_Sets(i,1) == Set_name(l))
tot_int_sets = tot_int_sets + 1;
if ObliqueToppling(i) == 1
Hc = Hc + 1;
critic_OT=critic_OT+1;
end
end
end
Set_OT(index_Set_OT,3) = tot_int_sets;
Set_OT(index_Set_OT,4) = Hc;
end
end
for i = 1 : numel(Set_OT(1,:))
if i ==1 || i==2
Set_OT(index_Set_OT+1,i)=NaN;
elseif i==3
Set_OT(index_Set_OT+1,i)=sum(Set_OT(:,i));
elseif i==4
Set_OT(index_Set_OT+1,i)=critic_OT;
end
end
clearvars OT_T
OT_T=table(Set_OT(:,1),Set_OT(:,2),Set_OT(:,3),Set_OT(:,4));
OT_T.Properties.VariableNames{'Var1'} = 'Set_i';
OT_T.Properties.VariableNames{'Var2'} = 'Set_j';
OT_T.Properties.VariableNames{'Var3'} = 'TotIntersections';
OT_T.Properties.VariableNames{'Var4'} = 'TotCritical';
%BasalPlane-Statistic
Set_BP=zeros(nSet+1,5);
critic_BP_Hc=zeros;
critic_BP_Lc=zeros;
for k = 1 : nSet
Hc=0;
Lc=0;
Set_BP(k,1) = Set_name(k);%set number value
Set_BP(k,2) = sum(Set(:) == (Set_name(k)));%total number of discontinuity of set
for i = 1 : num_disc%value(%) of critical discontinuity
if Set(i)== Set_name(k) && BasalPlane(i) == 2
Hc = Hc + 1;
critic_BP_Hc=critic_BP_Hc+1;
elseif Set(i)== Set_name(k) == 1 && (BasalPlane(i) == 1 ||BasalPlane(i) ==3)
Lc=Lc+1;
critic_BP_Lc=critic_BP_Lc+1;
end
end
Set_BP(k,4) =Hc ;
Set_BP(k,5) =Lc ;
Set_BP(k,3) =Set_BP(k,4) + Set_BP(k,3);
end
for i = 1 : numel(Set_BP(1,:))
if i == 1
Set_BP(nSet+1,i)=NaN;
elseif i==2
Set_BP(nSet+1,i)=sum(Set_BP(:,i));
elseif i==3
Set_BP(nSet+1,i)=(critic_BP_Hc+critic_BP_Lc);
elseif i==4
Set_BP(nSet+1,i)=critic_BP_Hc;
elseif i==5
Set_BP(nSet+1,i)=critic_BP_Lc;
end
end
clearvars BP_T
BP_T=table(Set_BP(:,1),Set_BP(:,2),Set_BP(:,3),Set_BP(:,4),Set_BP(:,5));
BP_T.Properties.VariableNames{'Var1'} = 'Set';
BP_T.Properties.VariableNames{'Var2'} = 'TotDiscontinuities';
BP_T.Properties.VariableNames{'Var3'} = 'TotBasalPlane';
BP_T.Properties.VariableNames{'Var4'} = 'SlidingBasalPlane';
BP_T.Properties.VariableNames{'Var5'} = 'NoSlidingBasalPlane';
% % %
disp('PlanarSlidingKynematicAnalysis')
disp(PS_T)
disp(' ')
disp('WedgeSlidingKynematicAnalysis')
disp(WS_T)
disp(' ')
disp('FlexuralTopplingKynematicAnalysis')
disp(FT_T)
disp(' ')
disp('DirectTopplingKynematicAnalysis')
disp(DT_T)
disp(' ')
disp('ObliqueTopplingKynematicAnalysis')
disp(OT_T)
disp(' ')
disp('BasalPlaneOfTopplingKynematicAnalysis')
disp(BP_T)
disp('CriticalPercentual')
disp('%PlanarSliding %WedgeSliding %FlexuralToppling %DirectToppling %ObliqueToppling')
disp([num2str(planeslope(1,1)),' ',num2str(planeslope(1,2)),' ',num2str(PS_T.TotCritical(end)/PS_T.TotDiscontinuities(end)),' ',...
num2str(WS_T.TotCritical(end)/WS_T.TotIntersections(end)),' ',...
num2str(FT_T.TotCritical(end)/FT_T.TotDiscontinuities(end)),' ',...
num2str(DT_T.TotCritical(end)/DT_T.TotIntersections(end)),' ',...
num2str(OT_T.TotCritical(end)/OT_T.TotIntersections(end))])
end
end