Deep, intermediate, and shallow earthquakes with magnitude M > 4.5
(for all data the depth of earthquake is known)

Modules

Weighted objective functions

Input: partition Π={(Subscript[π, 1]),⃛,(Subscript[π, k])} of the set 𝒜 ⊂ (ℝ^n), (in case of WFmin[] only the set 𝒜), centers (Subscript[c, 1]),⃛,(Subscript[c, k]) ∈ (ℝ^n);
Output: weighted objective function values: ℱ(Π), \[FormalCapitalF](c), 𝒢(Π);
Submodules: none needed

WF[PI_]:=Module[{m=Length[Flatten[PI,1]],c},
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
Return[Sum[PI[[j,s,1]]*Norm[PI[[j,s,2]]-c[[j]]]^2,{j,Length[PI]} ,{s,Length[PI[[j]]]}] ]
];

WFmin[A_,w_,z_]:=Module[{},
Sum[w[[i]]Min[Table[Norm[z[[j]]-A[[i]]]^2,{j,Length[z]}]],{i,Length[A]}]
]

WG[PI_]:=Module[{m=Length[Flatten[PI,1]],c,set,cc},
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}];
set=Flatten[PI,1];
cc=Mean[WeightedData[set[[All,2]],set[[All,1]]]];
Return[Sum[Length[PI[[j]]]*Norm[cc-c[[j]]]^2,{j,Length[PI]} ]]
];

Fcir[PI_,Circ_]:=Module[{},
Return[
Sum[Sum[Dist1[PI[[j,s]],Circ[[j,1]],Circ[[j,2]]],{s,Length[PI[[j]]]}]
,{j,Length[PI]}]
]]

Weighted k -  Means

WKMeans[A_,w_,z_]:=Module[{m=Length[A],k=Length[z],n=Length[A[[1]]],c0,Pod,imin,PI,It=0,tab},
Pod=Table[{w[[i]],A[[i]]},{i,m}]; 
c0=z; 
PI=Table[{},{i,k}];  
			(* Minimal distance principle *)
Do[
tab=Table[Norm[Pod[[i,2]]-c0[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
						(* New centers *)
While[
c=Table[Mean[WeightedData[PI[[j,All,2]],PI[[j,All,1]]]],{j,Length[PI]}]; 
		(* The case of reducing the number of clusters in a partition
				in the case n=2 - can be extended to an arbitrary n *)
If[Length[c]<Length[c0],
Print["There was a case of reducing the number of clusters in the partition"]];
Chop[Norm[c-c0]]≠0,c0=c;It=It+1;
				(* New partition *)
PI=Table[{},{i,k}]; 
Do[
tab=Table[Norm[Pod[[i,2]]-c0[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];  
               ];
{PI,c0,WF[PI],It}]

Distance between the point T∈(ℝ^n)and the circle Circ(c, r)in (ℝ^n)

Input: point T∈(ℝ^n), center c ∈ (ℝ^n)and radius r >0;
Output: Distance from the point T to the circle Circ(c, r);
Submodules: None
					(* LAD distance *)
Dist1[T_,c_,r_]:=Module[{},Return[Abs[Norm[T-c]-r] ]
                                             ];       
				(* TLS quasi distance *)
Dist2[T_,c_,r_]:=Module[{},Return[(Norm[T-c]-r)^2 ]
                                             ]; 
			(* Algebraic fit criterion  *)
DistA[T_,c_,r_]:=Module[{},Return[(Norm[T-c]^2-r^2)^2]
                                             ];

Sphere Representative

(* LAD(approach=1),  TLS(approach=2); Algebraic (approach=3); eps = accoracy *)
SphRep[A_,C0_,r0_,eps_,approach_]:=Module[{p1,q1,t1,c1,r1,F1,T,min,out},
Which[approach==1,
F[rr_]:=Sum[Dist1[A[[i]],C0,rr],{i,Length[A]}];
T=Timing[min=FindMinimum[F[rr],{rr,r0},Method->"Newton"]];
r1=rr/.min[[2]]; F1=T[[2,1]];
out={approach,r1,F1,T[[1]]},
approach==2,
F[rr_]:=Sum[Dist2[A[[i]],C0,rr],{i,Length[A]}];
T=Timing[min=FindMinimum[F[rr],{rr,r0},Method->"Newton"]];
r1=rr/.min[[2]]; F1=T[[2,1]];
out={approach,r1,F1,T[[1]]},
approach==3,
F[rr_]:=Sum[DistA[A[[i]],C0,rr],{i,Length[A]}];
T=Timing[min=FindMinimum[F[rr],{rr,r0},Method->"Newton"]];
r1=rr/.min[[2]];F1=T[[2,1]];
out={approach,r1,F1,T[[1]]}
];
out
]

KCS (k - means for spheres) - general

(* Ako je Ind≠0, daju se međukoraci sa slikama *)
KCS[A_,Circ0_,Iter_,Ind_]:=Module[{it=0,eps=0.05,m=Length[A],k=Length[Circ0],Sph,tab,p0,q0,t0,r0,fc,min,imin,F00,F0,F1},
Sph=Table[0,{j,k}];
PI=Table[{},{i,k}];
Do[
tab=Table[Dist1[A[[i]],Circ0[[j,1]],Circ0[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
F0=F00=Fcir[PI,Circ0]; 
Print["hat F=",F0, "; k= ",k, Table[Length[PI[[j]]],{j,k}]];
				 (* Računanje kružnica-centara LAD distance*)
Do[
r0=Median[Table[Norm[ C0-PI[[j,s]]],{s,Length[PI[[j]]]}]];
fc[rr_]:=Sum[Dist2[PI[[j,s]],C0,rr],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[rr],{rr,r0}];
Sph[[j]]={C0,rr/.min[[2]]};
,{j,k}];
F1=Fcir[PI,Sph] ;  
Print["Sph: ",Sph,"; F1 = ",F1];
While[
(F0-F1)/F00>eps||it<Iter,
it=it+1;
If[Ind≠0,
Print["It ",it,": Sphere: ", Sph,";  F = ",{F0,F1},"; del F = ",100(F0-F1)/F00];
sl5= Graphics3D[{{Opacity[.2],Sphere[Sph[[1,1]],Sph[[1,2]]]},{Opacity[.3],Sphere[Sph[[2,1]],Sph[[2,2]]]},{Opacity[.8],Sphere[Sph[[3,1]],Sph[[3,2]]]}}];
Print[Show[
ListPointPlot3D[PI,PlotStyle->{{Blue,PointSize->.013},{Red,PointSize->.013},{Green,PointSize->.013}}],sl5,slSph1,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]];];
F0=F1;
PI=Table[{},{j,k}];
Do[
tab=Table[Dist1[A[[i]],Sph[[j,1]],Sph[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
Sph=Table[0,{j,k}];
Do[
r0=Median[Table[Norm[ C0-PI[[j,s]]],{s,Length[PI[[j]]]}]];
fc[rr_]:=Sum[Dist1[PI[[j,s]],C0,rr],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[rr],{rr,r0}];
Sph[[j]]={C0,rr/.min[[2]]};
,{j,k}];
F1=Fcir[PI,Sph]; 
];
{PI,Sph,F1,it}
]

KCSO (k - closest concentric spheres with center O)

(* If Ind≠0, shows intermediate steps with pictures *)
KCSO[A_,rr0_,Iter_,Ind_]:=Module[{it=0,eps=0.00005,m=Length[A],k=Length[rr0],O={0,0,0},tab,r0,fc,min,imin,F00,F0,F1},
Sph=Table[0,{j,k}];
PI=Table[{},{i,k}];
Do[
tab=Table[Dist1[A[[i]],O,rr0[[j]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI]; 
F0=F00=Fcir[PI,Table[{O,rr0[[j]]},{j,k}]]; 
Print["hat F=",F0, "; k= ",k, Table[Length[PI[[j]]],{j,k}]];
				 (* Sphere centers - LAD distance *)
Do[
r0=Median[Table[Norm[ PI[[j,s]]],{s,Length[PI[[j]]]}]];
fc[rr_]:=Sum[Dist1[PI[[j,s]],O,rr],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[rr],{rr,r0}];
Sph[[j]]={O,rr/.min[[2]]};
,{j,k}];
F1=Fcir[PI,Sph] ;  
Print["Sph: ",Sph,"; F1 = ",F1];
While[
delF=(F0-F1)/F00;
delF>eps&&it<Iter,
it=it+1;
If[Ind≠0,
Print["It ",it,": Sphere: ", Sph,";  F = ",{F0,F1},"; delF = ",delF];
sl5= Graphics3D[{{Opacity[.2],Sphere[Sph[[1,1]],Sph[[1,2]]]},{Opacity[.3],Sphere[Sph[[2,1]],Sph[[2,2]]]},{Opacity[.8],Sphere[Sph[[3,1]],Sph[[3,2]]]}}];
Print[Show[
ListPointPlot3D[PI,PlotStyle->{{Blue,PointSize->.013},{Red,PointSize->.013},{Green,PointSize->.013}}],sl5,slSph1,PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]];];
F0=F1;
PI=Table[{},{j,k}];
Do[
tab=Table[Dist1[A[[i]],Sph[[j,1]],Sph[[j,2]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],A[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
Sph=Table[0,{j,k}];
Do[
r0=Median[Table[Norm[ PI[[j,s]]],{s,Length[PI[[j]]]}]];
fc[rr_]:=Sum[Dist1[PI[[j,s]],O,rr],{s,Length[PI[[j]]]}];
min=FindMinimum[fc[rr],{rr,r0}];
Sph[[j]]={O,rr/.min[[2]]};
,{j,k}];
F1=Fcir[PI,Sph]; 
];
{PI,Sph[[All,2]],F1,it}
]

Structure of strong earthquakes

(* Mag[M_,NM_,D1_,D2_]: M-magnituda (4 za 4≤M<5 ili 5 za M≥5), NM-broj potresa magnitude M, D1-donji rub sloja, D2-gornji rub sloja *)
Mag[M_,NM_,D1_,D2_]:=Module[{Nr=0},
Do[
If[D1≥  asd[[i,5]]≥D2&&asd[[i,6]]≥M,Nr=Nr+1]
,{i,Length[asd]}];
Print["Broj potresa magnitude ",M," u sloju [",D1,",",D2,"] = ",Nr," (",Round[100 Nr/NM//N,-.5],"%)"]
]

SymDIRECT

Modul for solving the global optimization problem for a simmetric Lipschitz continuous function. 
Mathematica codes written by I.Vazler;
Input: ordinary or compiled function, Domain ,
(*{Eps->10^(-4),MaxIter->20,DMin->10.^(-2)} *);
Output: optimal variables values, optimal objective function value , CPU-time, ;
Submodules: none needed

(* For a given center and length of the side of the hyperrectangle the function test returns:
  1 If the whole hyperrectangle is located in the region \Delta
0 If a hyperrectangle lies in the region \Delta only partially
-1 If a hyperrectangle appears outside the region \Delta (or if a hyperrectangle has only one point in \Delta)*)
test[c_,dd_]:=Module[{n=Length[c],S,d=dd/2,unutra=1},
S=Table[1,{n}];
Do[
If[c[[i]]-d[[i]]<c[[i+1]]+d[[i+1]],
unutra=0;Break[];
]
,{i,1,n-1}];
If[unutra==0,
Do[
If[c[[i]]+S[[i]]*d[[i]]<c[[i+1]]+d[[i+1]],
If[c[[i]]+S[[i]]*d[[i]]>c[[i+1]]-d[[i+1]],
S[[i+1]]=-1,
unutra=-1;Break[];
]
]
,{i,1,n-1}];
];
Return[unutra];
]

Options[SymDIRECT]={Eps->0,MaxIter->20,DMin->10.^(-2)};
SymDIRECT[f_,dom_,OptionsPattern[]]:=Module[{n,ff,c,d,dmax,fc,dmaxrazni,pot0,pot1,pot2,temp,temp2,ind,novic,novid,novifc,noviunutra,red,TMin,dtest,anim={},unutra,brisi},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
SymDIRECT::maxiter="The maximal number `1` of iterations has been reached.";
iter=0; Nf=0;
n=Length[dom];
(* scaling the function on the hypercube with length of the side 1 *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* initial center and diameters *)
c={Table[1/2,{i,n}]};
d={Table[1,{i,n}]};
unutra={test[c[[1]],d[[1]]]};
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@c[[1]]};   Nf=Nf+1;
TMin={c[[1]],fc[[1]]};
(* MAIN LOOP *)
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* STEP 1: searching for potentially optimal hyperrectangles *)
dmaxrazni=Union[dmax];
(* All data grouped according to d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* Data with smallest f grouped according to d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* Checking other conditions *)
pot2={pot1[[-1]]};
Do[temp=Min[Table[(fc[[pot1[[i,1]]]]-fc[[pot2[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot2[[j,1]]]]),{j,Length[pot2]}]];
If[temp>OptionValue[Eps]&&(i==1||Max[Table[(fc[[pot1[[i,1]]]]-fc[[pot1[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot1[[j,1]]]]),{j,i-1}]]<temp)
,AppendTo[pot2,pot1[[i]]]];
,{i,Length[dmaxrazni]-1,1,-1}];
pot2=Flatten[pot2,1];
(* STEP 2: dividing potentially optimal hyperrectangles *)
brisi={};
Do[
(* Dividing directions *)
ind=Flatten[Position[d[[pot]],Max[d[[pot]]]]];
novic=Table[c[[pot]],{2*Length[ind]}];
novid=Table[d[[pot]],{2*Length[ind]}];
d[[pot,ind]]/=3;
dmax[[pot]]/=3;
Do[
novic[[2i-1,ind[[i]]]]-=d[[pot,ind[[i]]]];
novic[[2i,ind[[i]]]]+=d[[pot,ind[[i]]]];
,{i,Length[ind]}];
novifc=Table[ff@@novic[[i]],{i,Length[novic]}];
Nf=Nf+Length[novic];
red=Ordering[Table[Min[{novifc[[2i-1]],novifc[[2i]]}],{i,Length[ind]}]];
Do[novid[[2red[[j]]-1,ind[[red[[i]]]]]]=novid[[2red[[j]],ind[[red[[i]]]]]]/=3,{i,1,Length[red]},{j,i,Length[red]}];
(* Eliminating those that are not inside *)
If[unutra[[pot]]==0,
(* Central hyperrectangle can go out or stay completely inside *)
If[test[c[[pot]],d[[pot]]]==-1,AppendTo[brisi,{pot}],unutra[[pot]]=test[c[[pot]],d[[pot]]]];
temp={};noviunutra={};
Do[
temp2=test[novic[[i]],novid[[i]]];
If[temp2==-1,AppendTo[temp,{i}],AppendTo[noviunutra,temp2]];
,{i,Length[novic]}];
novic=Delete[novic,temp];
novid=Delete[novid,temp];
novifc=Delete[novifc,temp];
,
(* If the divided hyperrectangle lies completely inside, all subhyperrectangles remain inside too *)
noviunutra=Table[1,{Length[novic]}];
];
c=Join[c,novic];
d=Join[d,novid];
dmax=Join[dmax,Map[Max,novid]];
fc=Join[fc,novifc];
unutra=Join[unutra,noviunutra];
(* Update of actual minimum *)
temp=Position[novifc,Min[novifc]][[1,1]];
If[novifc[[temp]]<TMin[[2]],TMin=N[{novic[[temp]],novifc[[temp]]},10]];
,{pot,pot2}];
(* Deletion of those hyperrectangles which went outside *)
c=Delete[c,brisi];
d=Delete[d,brisi];
dmax=Delete[dmax,brisi];
fc=Delete[fc,brisi];
unutra=Delete[unutra,brisi];
(*Print[iter,": ",{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];*)
If[iter==OptionValue[MaxIter],Message[SymDIRECT::maxiter,OptionValue[MaxIter]];];
(* Figure in 2D case *)
(*If[n==2 && iter<100,
AppendTo[anim,
Graphics[{{Green,Rectangle[]},Table[{Blue,Tooltip[Rectangle[c[[i]]-d[[i]]/2+d[[i]]/10,c[[i]]+d[[i]]/2-d[[i]]/10],unutra[[i]]]},{i,Length[d]}],Point[c]}]
];
];*)
];
,{iter,Length[pot2]}];
(* Back to the initial domain *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* Images printing in 2D case *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

Data + Initial approximation

(* 4.5 < M < 5 *)  
SetDirectory[NotebookDirectory[]]; 
asd=Import["pod4.txt","Table"];
total45=Length[asd];
priv={};
Do[
If[asd[[i,6]]≠ 5&&asd[[i,5]]≠0,priv=Append[priv,asd[[i]]]];
,{i,Length[asd]}];
asd=priv;
phi45=asd[[All,3]]; 
th45=asd[[All,4]];  
depth45=-asd[[All,5]]; 
Print["EQ 4-5:  TotalEQ = ",total45,"; Total (ne) = ", Length[depth45]]

(* 5 ≤ M < 6 *)
asd=Import["5-6.txt","Table"];
total56=Length[asd];
priv={};
Do[
If[asd[[i,6]]≠ 6&&asd[[i,5]]≠0,priv=Append[priv,asd[[i]]]];
,{i,Length[asd]}];
asd=priv;
phi56=asd[[All,3]]; 
th56=asd[[All,4]];  
depth56=-asd[[All,5]]; 
Print["EQ 5-6:  TotalEQ = ",total56,"; Total (ne) = ", Length[depth56]]

(* 6 ≤ M < 7 *)
asd=Import["6-7.txt","Table"];
total67=Length[asd];
priv={};
Do[
If[asd[[i,6]]≠7&&asd[[i,5]]≠0,priv=Append[priv,asd[[i]]]];
,{i,Length[asd]}];
asd=priv;
phi67=asd[[All,3]]; 
th67=asd[[All,4]];  
depth67=-asd[[All,5]]; 
Print["EQ 6-7:  TotalEQ = ",total67,"; Total (ne) = ", Length[depth67]]

(* 7 ≤ M  *)
asd=Import["7-.txt","Table"];
total7=Length[asd];
priv={};
Do[
If[asd[[i,6]]≠7&&asd[[i,5]]≠0,priv=Append[priv,asd[[i]]]];
,{i,Length[asd]}];
asd=priv;
phi7=asd[[All,3]]; 
th7=asd[[All,4]];  
depth7=-asd[[All,5]]; 
Print["EQ 7:  TotalEQ = ",total7,"; Total (ne) = ", Length[depth7]]
Print["All EQ: ",sve={total45,total56,total67,total7}," = ",Total[sve]]
Print["EQ without 0: ",sve={Length[phi45],Length[phi56],Length[phi67],Length[phi7]}," = ",Total[sve]]

205886-186993

Transformation on [0, 1]

C0={0,0,0}; Print["C0 = ",C0];
phi=Join[phi45,phi56,phi67,phi7]; phi=Pi phi/180.;
th=Join[th45,th56,th67,th7]; th=Pi th/180.;
depth=Join[depth45,depth56,depth67,depth7];  
m=Length[depth]; Print["m = ",m];
 min=Min[depth]; (* najdublje *)
 max=Max[depth]; (* najpliće *)
Print["{najdublje, najpliće}: ",{min, max}," -> [0,1]"]
lin[t_]:=-t/min+1.
lin1[t_]=InverseFunction[lin][t];
B=lin[depth];
II={0,1};
w=Table[1,{i,Length[depth]}];
n=1;k=3;
line=Graphics[{Thickness[0.004], Line[{{0,0},{1,0}}]}];
slB=ListPlot[Table[{B[[j]],0},{j,m}],PlotStyle->{PointSize[.01]}];
slB=Show[line,slB,PlotRange->{II,{-.05,.05}},Axes->{True,False},Ticks->{None,None},AspectRatio->Automatic,ImageSize->Medium]
						(* The set A *)
A=Table[B[[i]]{Cos[phi[[i]]] Cos[th[[i]]],Cos[phi[[i]]] Sin[th[[i]]],Sin[phi[[i]]]},{i,m}];
slA=ListPointPlot3D[A,PlotStyle->{PointSize->.01} ];
Print[Style["{Norm[A[[j]]] = B[[j]]} for all j !!!",Blue,16]];

Are there outliers in B?

dubine=-depth ;
med=Median[dubine];
mad=1.483Median[Abs[dubine-med]]; Print["MAD=",mad];
dubinetil=Abs[dubine-med]/mad;
Outliers={};
Do[
If[dubinetil[[i]]>2.5,Outliers=Append[Outliers,{i,dubine[[i]]}]]
,{i,Length[dubinetil]}]
Outliers
Print["Od ukupno ",m," podataka postoji ",Length[Outliers]," outliers, tj., ",100.Length[Outliers]/m,"% outliers"]

Example: outliers determining

pod={9.05,2.83,3.00,3.16,4.12,3.00,3.50};
med=Median[pod]
Abs[pod-med]
mad=1.483Median[Abs[pod-med]]; Print["MAD=",mad];
podtil=Abs[pod-med]/mad

One Sphere

Print[Style["Approach 1:  r = ",Blue,14],r1=Median[B],"; Actual depth = ", lin1[r1],"km"]
Print[Style["Approach 2:  r = ",Blue,14],r2=Mean[B],"; Actual depth = ", lin1[r2]," km"]
Print[Style["Approach A:  r = ",Blue,14],rA=Sqrt[Mean[Table[B[[s]]^2,{s,m}]]],"; Actual depth = ", lin1[rA]," km"]

Figure

r0=r1; Print["r0 = ", r0]; 
sph=SphRep[A,C0,r0,-.0005,1];  (* Birao sam LAD distance "1" *)
Print[Style["Approach: ",Blue,14],sph[[1]],"\nr = ",sph[[2]],";  F=",sph[[3]],"; CPU=",sph[[4]]]
AG=Table[{-1,1},3];
slSph1=Graphics3D[{LightBlue,Opacity[.1],Sphere[C0,1]}];
slSph=Graphics3D[{LightBlue,Opacity[.3],Sphere[C0,sph[[2]]]}];
slc=ListPointPlot3D[{C0}, PlotStyle->{Brown,PointSize[.06],Opacity[.8]}];
Show[slA,slSph,slSph1,
PlotRange->AG,BoxRatios->{1,1,1},AxesLabel->{x,y,z},ImageSize->Small]
slcSym=ListPlot[{{r0,0}},PlotMarkers->{"✶",15},PlotStyle->{Red,PointSize[.02],Opacity[.8]}(*,PlotRange->{II,{-.05,.05}},Axes->{True,False},Ticks->{{{0,"0"},{cenSymDir[[1]],"c1"},{cenSymDir[[2]],"c2"},{cenSymDir[[3]],"c3"},{1,"1"}},None},AspectRatio->Automatic*)];
Show[slB,slcSym,ImageSize->300]
Print["Stvarna dubina = ",lin1[sph[[2]]]];

Three Sphere

Initial approximation for k=3

k=3;
cFunk0[var_]:=Compile[var,Evaluate[Sum[Min[Table[Norm[Table[(Subscript[x, j]),{j,1+(jk-1)n,1+jk n-1}]-B[[i]]]^2,{jk, 1,k}]],{i,m}]]];
FFCompile=cFunk0[Table[{(Subscript[x, i]),_Real},{i,n k}]];

(* SymDIRECT *)
Print[Style["SymDIRECT:",16,Blue]];
dom=Table[II,{k}];
tSym=Timing[
minSym=SymDIRECT[FFCompile,dom,DMin->10^(-4),MaxIter->30,Eps->10^(-4)]
];
cenSymDir=minSym[[1]];
Print[Style["SymDIRECT Centers: ",Green,16],cenSymDir];
Print["F= ",minSym[[2]],"; IT = ",iter,"; CPU = ",tSym[[1]],"; Nf = ",Nf];
slcSym=ListPlot[Table[{cenSymDir[[j]],0},{j, k}],PlotMarkers->{"✶",15},PlotStyle->{Red,PointSize[.02],Opacity[.8]}(*,PlotRange->{II,{-.05,.05}},Axes->{True,False},Ticks->{{{0,"0"},{cenSymDir[[1]],"c1"},{cenSymDir[[2]],"c2"},{cenSymDir[[3]],"c3"},{1,"1"}},None},AspectRatio->Automatic*)];
Show[slB,slcSym,ImageSize->300]

AG=Table[{-1.,1.},3];
Circ0=Table[{C0,cenSymDir[[j]]},{j,3}]
kcs=Timing[KCSO[A,cenSymDir,2,1]];
Print["r: ",rd=kcs[[2,2]]];
Print["Actual Depths: ",Rd=-lin1[kcs[[2,2]]] ];
Print["CPU(KCS) = ",kcs[[1]]," seconds = ", kcs[[1]]/60," min"]

Analysis of results in the Unit Sphere model

Print[Style["Radii of spheres inside Unit Sphere: ",Blue,14],rd];
ed={lin[max],Mean[{rd[[1]],rd[[2]]}],Mean[{rd[[2]],rd[[3]]}],lin[min]};
Print[Style["Radii of boundary spheres",Blue,14],Style["(they delimit clusters of transformed depths within the Unit Sphere):\n",Gray,12],ed];
Print[Style["Spherical regions inside Unit Sphere with the number of earthquake events:",Blue,14]];Print["Deep: ",{ed[[1]],ed[[2]]}," | ",Length[kcs[[2,1,1]]]];
j=1;list=Table[Norm[kcs[[2,1,j,s]]],{s,Length[kcs[[2,1,j]]]}];
Print["Mean[π(",j,")] = ",Mean[list],"; StandardDeviation = ",StandardDeviation[list]]
Print["Intermediate: ",{ed[[2]],ed[[3]]}," | ",Length[kcs[[2,1,2]]]];
j=2;list=Table[Norm[kcs[[2,1,j,s]]],{s,Length[kcs[[2,1,j]]]}];
Print["Mean[π(",j,")] = ",Mean[list],"; StandardDeviation = ",StandardDeviation[list]]
Print["Shallow: ",{ed[[3]],ed[[4]]}," | ",Length[kcs[[2,1,3]]]]
j=3;list=Table[Norm[kcs[[2,1,j,s]]],{s,Length[kcs[[2,1,j]]]}];
Print["Mean[π(",j,")] = ",Mean[list],"; StandardDeviation = ",StandardDeviation[list]]

Analysis of real results

ed1=-lin1[ed];
Print[Style["Depths (in kilometers) of cluster centers in the Geoid: ",Blue,14],Round[-lin1[rd],-.005]]
Print[Style["Radii of boundary spheres (in km): ",Blue,14],ed1]
Print[Style["Depth ranges (in kilometers) with numbers of earthquake events: ",Blue,14]];Print["Shallow:",{ ed1[[1]], ed1[[2]]}," | ",Length[kcs[[2,1,1]]]];
j=1;list=Table[-lin1[Norm[kcs[[2,1,j,s]]]],{s,Length[kcs[[2,1,j]]]}];
Print["Mean[π(",j,")] = ",Mean[list],"; StandardDeviation = ",StandardDeviation[list]]
Print["Intermediate: ",{ ed1[[2]], ed1[[3]]}," | ",Length[kcs[[2,1,2]]]]
j=2;list=Table[-lin1[Norm[kcs[[2,1,j,s]]]],{s,Length[kcs[[2,1,j]]]}];
Print["Mean[π(",j,")] = ",Mean[list],"; StandardDeviation = ",StandardDeviation[list]]
Print["Deep: ",{ ed1[[3]], ed1[[4]]}," | ",Length[kcs[[2,1,3]]]]
j=3;list=Table[-lin1[Norm[kcs[[2,1,j,s]]]],{s,Length[kcs[[2,1,j]]]}];
Print["Mean[π(",j,")] = ",Mean[list],"; StandardDeviation = ",StandardDeviation[list]]

Structure of strong earthquakes

Mag[5,72137,24,1]
Mag[5,72137,348.5,24]
Mag[5,72137,663,348.5]

Mag[6,5296,24,1]
Mag[6,5296,348.5,24]
Mag[6,5296,663,348.5]