load ('skin.mat');
skin_tri=double(currytri');
skin_pnt=double(curryloc');

load ('innerSkull.mat');
innerSkull_tri=double(currytri');
innerSkull_pnt=double(curryloc');

load ('outerSkull.mat');
outerSkull_tri=double(currytri');
outerSkull_pnt=double(curryloc');


tr1=TriRep(skin_tri,skin_pnt(:,1),skin_pnt(:,2),skin_pnt(:,3));
tr2=TriRep(outerSkull_tri,outerSkull_pnt(:,1),outerSkull_pnt(:,2),outerSkull_pnt(:,3));
tr3=TriRep(innerSkull_tri,innerSkull_pnt(:,1),innerSkull_pnt(:,2),innerSkull_pnt(:,3));
     
figure (1)

       trisurf(tr1,'facealpha',0.1,'facecolor',[0 0 1]); hold on
       trisurf(tr2,'facealpha',0.2,'facecolor',[0 1 0]); hold on
       trisurf(tr3,'facecolor',[1 0 0])
        
        %% BEM Computation-calculating the related electromagnetic model

Cskin_surface= 0.33; %scalp
Couter_skull_surface=0.0042; %skull
Cinner_skull_surface=0.33; %brain
vol=[];

vol.bnd(1).tri=skin_tri;
vol.bnd(1).pnt=skin_pnt;

vol.bnd(2).tri=outerSkull_tri;
vol.bnd(2).pnt=outerSkull_pnt;

vol.bnd(3).tri=innerSkull_tri;
vol.bnd(3).pnt=innerSkull_pnt;

vol.cond=[Cskin_surface Couter_skull_surface Cinner_skull_surface];

cfg=[];
cfg.method         = 'openmeeg';
vol = ft_prepare_headmodel(cfg, vol);