% function [E_X,E_Y,E_Z]=E_FIELD _CALC(POS)
% POS
clear
clc

POS_Fix=zeros(4,3);
% point3(0,5,10.5),point3(5,0,11.83),point3(0,-5,13.17),point3(-5,0,14.5)
X_Fix=[0 5 0 -5]; Y_Fix=[5 0 -5 0]; Z_Fix=[10.5 11.83 13.17 14.5];
for i=1:4
POS_Fix(i,1)=X_Fix(i);
POS_Fix(i,2)=Y_Fix(i);
POS_Fix(i,3)=Z_Fix(i);
end

% X_Var=[0 0 0 0 0]; Y_Var=[0 0 0 0 0]; Z_Var=[40 20 0 -20 -40];
% X_Var=[0 -25 -25 -25 20 25 0 0 0 0 10 9 10 15 -5 -10 -15 -20 -25 0 0 -25 -25 -25 20 25 0 0 0 0 10 9 10 15 -5 -10 -15 -20 -25 0 ];
% Y_Var=[0 0 0 0 0 0 0 0 10 20 20 -20 -23 -15 -5 -13 -20 -10 -4 -20 0 0 0 0 0 0 0 0 10 20 20 -20 -23 -15 -5 -13 -20 -10 -4 -20 ];
% Z_Var=[40 50 30 35 35 30 30 35 40 33 45 45 40 45 30 27 30 35 40 35 -40 -50 -30 -35 -35 -30 -30 -35 -40 -33 -45 -45 -40 -45 -30 -27 -30 -35 -40 -35 ];

% X_Var(24)=25;

X_Var=[12,15,11,11,-12,-11,-7,6,-12,18,-24,25,0,-9,-6,20,-19,21,-14,12,-26,19,3,5,-4,-17,-15,16,19,1,22,-7,-2,7,18,4,-1,-14,2,13];
Y_Var=[18,10,-12,-10,-8,-6,-3,-10,0,6,11,9,20,-15,-4,-10,15,16,-2,7,10,18,-27,-9,-10,-16,-13,-2,6,15,-12,27,16,5,-8,-13,12,24,11,9];
Z_Var=[28,-49,-42,-37,-33,-29,-31,51,33,38,32,-43,38,-35,-41,52,-26,-34,44,-31,-40,-44,36,-42,36,-38,32,32,-50,48,49,29,34,27,-37,-44,41,-48,-44,-33];
    
POS_Var=zeros(40,3);
for i=1:40
POS_Var(i,1)=X_Var(i);
POS_Var(i,2)=Y_Var(i);
POS_Var(i,3)=Z_Var(i);
end

% POS_Var=POS;
% POS_Var=zeros(10,3);
% for i=1:10
% POS_Var(i,2)=15-3*i;
% end
% POS_Var(:,1)=50;
% POS_Var=POS_Var*1e0;

clear fem
fem.const={'Q1', 1.3e-19};
clear draw
% draw.s.objs={cylinder3(30,33,'pos',[0 0 25])+cylinder3(30,33,'pos',[0 0 -58])+cylinder3(4,50,'pos',[0 0 -25])+cylinder3(5.5,50,'pos',[0 0 -25])};
% draw.s.name={'C01'};
mm=0; nn=5; pp=10.4;
 g_fix=sphere3(0.1e0,'pos',[POS_Fix(1,1) POS_Fix(1,2) POS_Fix(1,3)]);%,'axis',[0 0 1],'rot',0);
for i=2:4
g1=sphere3(0.1e0,'pos',[POS_Fix(i,1) POS_Fix(i,2) POS_Fix(i,3)]);%,'axis',[0 0 1],'rot',0);
g_fix=geomcomp(g_fix,g1);
end

R=0.99;

 g_var1=sphere3(R,'pos',[POS_Var(1,1) POS_Var(1,2) POS_Var(1,3)]);%,'axis',[0 0 1],'rot',0);
for i=2:10
g_1=sphere3(R,'pos',[POS_Var(i,1) POS_Var(i,2) POS_Var(i,3)]);%,'axis',[0 0 1],'rot',0);
% g_var1=g_var1+g_ 1;
g_var1=geomcomp(g_var1,g_1);
end
 g_var2=sphere3(R,'pos',[POS_Var(11,1) POS_Var(11,2) POS_Var(11,3)]);%,'axis',[0 0 1],'rot',0);
for i=12:20
g_2=sphere3(R,'pos',[POS_Var(i,1) POS_Var(i,2) POS_Var(i,3)]);%,'axis',[0 0 1],'rot',0);
g_var2=geomcomp(g_var2,g_2);
end
% g_var2={sphere3(R,'pos',[POS_Var(11,1) POS_Var(11,2) POS_Var(11,3)],'axis',[0 0 1],'rot',0)+...
%     sphere3(R,'pos',[POS_Var(12,1) POS_Var(12,2) POS_Var(12,3)],'axis',[0 0 1],'rot',0)+...
%     sphere3(R,'pos',[POS_Var(13,1) POS_Var(13,2) POS_Var(13,3)],'axis',[0 0 1],'rot',0)+...
%     sphere3(R,'pos',[POS_Var(14,1) POS_Var(14,2) POS_Var(14,3)],'axis',[0 0 1],'rot',0)+...
%     sphere3(R,'pos',[POS_Var(15,1) POS_Var(15,2) POS_Var(15,3)],'axis',[0 0 1],'rot',0)+...
%     sphere3(R,'pos',[POS_Var(16,1) POS_Var(16,2) POS_Var(16,3)],'axis',[0 0 1],'rot',0)+...
%     sphere3(R,'pos',[POS_Var(17,1) POS_Var(17,2) POS_Var(17,3)],'axis',[0 0 1],'rot',0)+...
%     sphere3(R,'pos',[POS_Var(18,1) POS_Var(18,2) POS_Var(18,3)],'axis',[0 0 1],'rot',0)+...
%     sphere3(R,'pos',[POS_Var(19,1) POS_Var(19,2) POS_Var(19,3)],'axis',[0 0 1],'rot',0)+...
%     sphere3(R,'pos',[POS_Var(20,1) POS_Var(20,2) POS_Var(20,3)],'axis',[0 0 1],'rot',0)};
    

 g_var3=sphere3(R,'pos',[POS_Var(21,1) POS_Var(21,2) POS_Var(21,3)]);%,'axis',[0 0 1],'rot',0);
for i=22:30
%     i
g_3=sphere3(R,'pos',[POS_Var(i,1) POS_Var(i,2) POS_Var(i,3)]);%,'axis',[0 0 1],'rot',0);
g_var3=geomcomp(g_var3,g_3);
end
 g_var4=sphere3(R,'pos',[POS_Var(31,1) POS_Var(31,2) POS_Var(31,3)]);%,'axis',[0 0 1],'rot',0);
for i=32:40
%     i
g_4=sphere3(R,'pos',[POS_Var(i,1) POS_Var(i,2) POS_Var(i,3)]);%,'axis',[0 0 1],'rot',0);
g_var4=geomcomp(g_var4,g_4);
end
% g_cyl=cylinder3(30,33,'pos',[0 0 25])+cylinder3(30,33,'pos',[0 0 -58])+cylinder3(4,50,'pos',[0 0 -25])+cylinder3(5.5,50,'pos',[0 0 -25]);

% g_cyl _ud=cylinder3(30,33,'pos',[0 0 -58])+cylinder3(30,33,'pos',[0 0 25]);
g_cyl_ud=cylinder3(40,33,'pos',[0 0 -58])+cylinder3(40,33,'pos',[0 0 25]);
g_cyl_c1=cylinder3(5.5,50,'pos',[0 0 -25]);
g_cyl_c2=cylinder3(4,50,'pos',[0 0 -25]);

% s.objs={cylinder3(30,33,'pos',[0 0 25])+cylinder3(30,33,'pos',[0 0 -58])+cylinder3(4,50,'pos',[0 0 -25])+cylinder3(5.5,50,'pos',[0 0 -25]),g_fix,g_var};

% s.objs={g_cyl _ud,g_cyl _c1,g_cyl _c2,g_fix,g_var1,g_var2,g_var3,g_var4};
% s.objs={g_cyl _ud,g_cyl _c2,g_cyl _c1,g_fix,g_var3};%,g_var2,g_var3,g_var4};%,g_var2};%,g_var3,g_var4};
s.objs={g_fix,g_var1,g_var2,g_var3,g_var4};%,g_var2,g_var3,g_var4};%,g_var2};%,g_var3,g_var4};

% s.objs={g_fix,g_var,g_cyl};
% s1.objs={g_cyl};
% s2.objs={g_fix};
% s3.objs={g_var};
% 

% s.objs={cylinder3(30e-10,33e-10,'pos',[0 0 25e-10])+cylinder3(30e-10,33e-10,'pos',[0 0 -58e-10])+cylinder3(4e-10,50e-10,'pos',[0 0 -25e-10])+cylinder3(5.5e-10,50e-10,'pos',[0 0 -25e-10]),g_fix,g_var};
% fem.draw=struct('s',s);


% draw.p.objs={point3(0,5,10.5),point3(5,0,11.83),point3(0,-5,13.17),point3(-5,0,14.5),point3(0,-5,10.5),point3(5,0,10.5),point3(-5,0,10.5),point3(0,5,11.83),point3(-5,0,11.83),point3(0,-5,11.83),point3(-5,0,13.17),point3(5,0,13.17),point3(0,5,13.17),point3(5,0,14.5),point3(0,5,14.5),point3(0,-5,14.5),point3(mm,nn,pp)};
% draw.p.name={'PT1','PT2','PT3','PT4','PT5','PT6','PT7','PT8','PT9','PT10','PT11','PT12','PT13','PT14','PT15','PT16','PT17'};
% p.objs={point3(0,5,10.5),point3(5,0,11.83),point3(0,-5,13.17),point3(-5,0,14.5),point3(0,-5,10.5),point3(5,0,10.5),point3(-5,0,10.5),point3(0,5,11.83),point3(-5,0,11.83),point3(0,-5,11.83),point3(-5,0,13.17),point3(5,0,13.17),point3(0,5,13.17),point3(5,0,14.5),point3(0,5,14.5),point3(0,-5,14.5),point3(mm,nn,pp)};
% fem.draw=struct('s',s,'p',p);

fem.draw=struct('s',s);



% fem.draw=draw;
% [g,st,stx]=geomcsg(fem);

% fem.geom=geomcsg(fem);
[fem.geom,st]=geomcsg(fem); 

% ind_var=find(st(:,5)==1);
% ind_fix=find(st(:,4)==1);
% ind_ud=find(st(:,1)==1 & st(:,5)==0);
% ind_div=find(st(:,2)==1 & st(:,5)==0 & st(:,4)==0);
% ind_cen=find(st(:,3)==1 & st(:,5)==0 & st(:,4)==0);
% 
% [mm,nn]=size(st);
% aei_ep=zeros(1,mm);
% aei_ep(ind_var')=1;
% aei_ep(ind_fix')=2;
% aei_ep(ind_ud')=3;
% aei_ep(ind_div')=4;
% aei_ep(ind_cen')=5;
% 
% % for i=1:mm
% %     find(find_fix==
% %     aei_ep(i)=
% 
% geomplot(fem, 'pointlabels', 'on');
% clear appl
% appl.dim= {'V'}
% appl.mode.class='EmElectrostatics';
% appl.name= 'emes';
% appl.equ.epsilonr={'1' '2' '3' '4' '5'};
% % appl.equ.ind = [1 1 1 1 2 2 2 2 3 3 3 3 3];
% appl.equ.ind = aei_ep;
% appl.equ.rho={'0' '10e-6' '20e-6' '30e-6' '40e-6'};
% % appl.equ.ind = [1 1 1 1 2 2 2 2 3 3 3 3 3];
% appl.equ.ind = aei_ep;
% 
% % appl.pnt.Q0={'0' 'Q1'};
% % appl.pnt.ind=[1 1 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 1 1 2 1];
% % appl.pnt.ind=[1 1 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 1 1 ];
% fem.appl=appl;
% fem=multiphysics(fem);
% fem.mesh=meshinit(fem);
% fem.xmesh=meshextend(fem);
% 
% 
% fem.sol=femlin(fem);
% % B=[-1:0.1:1;-1:0.1:1;-1:0.1:1];
% % E_X=postinterp(fem,'Ex',POS_Var');
% % E_Y=postinterp(fem,'Ey',POS_Var');
% % E_Z=postinterp(fem,'Ez',POS_Var');
% 
%