g1=line1([0.0075*sin(0.25*pi/180) -0.0105*sin(2.25*pi/180) -0.0105*sin(2.75*pi/180) -0.0075*sin(0.25*pi/180)], ...
[0.0075*cos(0.25*pi/180) 0.0105*cos(2.25*pi/180) 0.0105*cos(2.75*pi/180) 0.0075*cos(0.25*pi/180)]);
g2=arc1(0,0,0.0075,0.25*pi/180,(90-0.25)*pi/180);
g3=line1([-0.0075*cos(0.25*pi/180) -0.0105*cos(0.25*pi/180) -0.0105*cos(0.25*pi/180) -0.0075*cos(0.25*pi/180)], ...
[0.0075*sin(0.25*pi/180) 0.0105*sin(0.25*pi/180) -0.0105*sin(0.25*pi/180) -0.0075*sin(0.25*pi/180)]);
g4=arc1(0,0,0.0075,90.25*pi/180,(180-0.25)*pi/180);
g5=line1([-0.0075*sin(0.25*pi/180) -0.0105*sin(0.25*pi/180) 0.0105*sin(0.25*pi/180) 0.0075*sin(0.25*pi/180)], ...
[-0.0075*cos(0.25*pi/180) -0.0105*cos(0.25*pi/180) -0.0105*cos(0.25*pi/180) -0.0075*cos(0.25*pi/180)]);
g6=arc1(0,0,0.0075,180.25*pi/180,(270-0.25)*pi/180);
g7=line1([0.0075*cos(0.25*pi/180) 0.0105*cos(0.25*pi/180) 0.0105*cos(0.25*pi/180) 0.0075*cos(0.25*pi/180)], ...
[-0.0075*sin(0.25*pi/180) -0.0105*sin(0.25*pi/180) 0.0105*sin(0.25*pi/180) 0.0075*sin(0.25*pi/180)]);
g8=arc1(0,0,0.0075,270.25*pi/180,(360-0.25)*pi/180);
g9=geomcoerce('solid',{g1,g2,g3,g4,g5,g6,g7,g8});
g10=circ2(0.0115);
g11=circ2(0.0125);
g12=g11-g10;
Nsteps=200;
delta_phi=360/Nsteps;
for k=0:Nsteps
  flclear fem
  fem.const = {'freq','35.8[kHz]'};
  g13=rotate(g9,k*delta_phi*pi/180,0,0);
  g14=g10-g13;
  clear s
  s.objs={g12,g13,g14};
  s.name={'stator','rotor','air gap'};
  s.tags={'g12','g13','g14'};
  fem.draw=struct('s',s);
  fem.geom=geomcsg(fem);
  fem.mesh=meshinit(fem,'hauto',5,'methodsub','tri');
  fem.mesh=meshrefine(fem,'mcase',0,'rmethod','regular');
  fem.mesh=meshrefine(fem,'mcase',0,'rmethod','regular');
  clear appl
  appl.mode.class = 'AcoPlaneStrain';
  appl.module = 'ACO';
  appl.gporder = 4;
  appl.cporder = 2;
  appl.sshape = 2;
  appl.assignsuffix = '_acpn';
  clear prop
  prop.analysis='freq';
  appl.prop = prop;
  clear pnt
  pnt.Fy = {0,67};
  numb_of_ent=geominfo(fem.geom,'out','no','od',[0,1]);
  numb_of_pnt=numb_of_ent(1);
  numb_of_bnd=numb_of_ent(2);
  coord_of_pnt=geominfo(fem.geom,'out','mp');
  y_of_pnt=coord_of_pnt(2,:);
  x_of_pnt=coord_of_pnt(1,:);
  r_of_pnt=sqrt(x_of_pnt.^2+y_of_pnt.^2);
  load_pnt=find(y_of_pnt==0.0125);
  pnt.ind=ones(1,numb_of_pnt);
  pnt.ind(1,load_pnt)=2;
  appl.pnt = pnt;
  clear bnd
  bnd.name = {'','Fluid load'};
  bnd.Fx = {0,'nx_acpr*p'};
  bnd.Fy = {0,'ny_acpr*p'};
  ext_pnt1=load_pnt;
  ext_pnt2=find(y_of_pnt==-0.0125);
  int_pnt1=find(y_of_pnt==0.0115);
  int_pnt2=find(y_of_pnt==-0.0115);
  adj_matrix1=geominfo(fem.geom,'out','adj','odp',[0;1]);
  ext_bnd1=find(abs(sign(adj_matrix1{1}(ext_pnt1,:)))==1);
  ext_bnd2=find(abs(sign(adj_matrix1{1}(ext_pnt2,:)))==1);
  ext_bnd=[ext_bnd1,ext_bnd2];
  int_bnd1=find(abs(sign(adj_matrix1{1}(int_pnt1,:)))==1);
  int_bnd2=find(abs(sign(adj_matrix1{1}(int_pnt2,:)))==1);
  int_bnd=[int_bnd1,int_bnd2];
  bnd.ind=ones(1,numb_of_bnd);
  bnd.ind(1,int_bnd(1,1))=2;
  bnd.ind(1,int_bnd(1,2))=2;
  bnd.ind(1,int_bnd(1,3))=2;
  bnd.ind(1,int_bnd(1,4))=2;
  appl.bnd = bnd;
  clear equ
  equ.nu = {0.36,0.33};
  equ.thickness = {0.003,1};
  equ.rho = {8600,7850};
  equ.E = {9.8e10,2.0e11};
  equ.name = {'Solid domain','Fluid domain'};
  equ.usage = {1,0};
  adj_matrix2=geominfo(fem.geom,'out','adj','odp',[2;2]);
  for i=1:3
    l1(i)=length(find(abs(sign(adj_matrix2{1}(i+1,2:4)))==1));
  end
  fluid_pos=find(l1==2);
  equ.ind=ones(1,3);
  equ.ind(1,fluid_pos)=2;
  appl.equ = equ;
  appl.var = {'freq','freq'};
  fem.appl{1} = appl;
  clear appl
  appl.mode.class = 'AcoPressure';
  appl.module = 'ACO';
  appl.sshape = 2;
  appl.assignsuffix = '_acpr';
  clear prop
  clear weakconstr
  weakconstr.value = 'off';
  weakconstr.dim = {'lm3'};
  prop.weakconstr = weakconstr;
  appl.prop = prop;
  clear bnd
  bnd.type = {'SH','NA','P0'};
  bnd.name = {'','Structural acceleration',''};
  bnd.nacc = {0,'nx_acpn*u_tt_acpn+ny_acpn*v_tt_acpn',0};
  bnd.ind=ones(1,numb_of_bnd);
  bnd.ind(1,int_bnd(1,1))=2;
  bnd.ind(1,int_bnd(1,2))=2;
  bnd.ind(1,int_bnd(1,3))=2;
  bnd.ind(1,int_bnd(1,4))=2;
  bnd.ind(1,ext_bnd(1,1))=3;
  bnd.ind(1,ext_bnd(1,2))=3;
  bnd.ind(1,ext_bnd(1,3))=3;
  bnd.ind(1,ext_bnd(1,4))=3;
  appl.bnd = bnd;
  clear equ
  equ.name = {'Solid domain','Fluid domain'};
  equ.usage = {0,1};
  equ.ind=ones(1,3);
  equ.ind(1,fluid_pos)=2;
  appl.equ = equ;
  appl.var = {'freq','freq'};
  fem.appl{2} = appl;
  fem.frame = {'ref'};
  clear units;
  units.basesystem = 'SI';
  fem.units = units;
  clear ode
  clear units;
  units.basesystem = 'SI';
  ode.units = units;
  fem.ode=ode;
  fem=multiphysics(fem);
  fem.xmesh=meshextend(fem);
  fem.sol=femstatic(fem,'solcomp',{'v','u','p'},'outcomp',{'v','u','p'},'blocksize','auto','nonlin','off');
  blade_pnt=find(0.0104<=r_of_pnt<=0.0106);
  blade_bnd1=find(abs(sign(adj_matrix1{1}(blade_pnt(1),:)))==1);
  blade_bnd2=find(abs(sign(adj_matrix1{1}(blade_pnt(2),:)))==1);
  blade_bnd3=find(abs(sign(adj_matrix1{1}(blade_pnt(3),:)))==1);
  blade_bnd4=find(abs(sign(adj_matrix1{1}(blade_pnt(4),:)))==1);
  blade_bnd5=find(abs(sign(adj_matrix1{1}(blade_pnt(5),:)))==1);
  blade_bnd6=find(abs(sign(adj_matrix1{1}(blade_pnt(6),:)))==1);
  blade_bnd7=find(abs(sign(adj_matrix1{1}(blade_pnt(7),:)))==1);
  blade_bnd8=find(abs(sign(adj_matrix1{1}(blade_pnt(8),:)))==1);
  blade_bnd=[blade_bnd1,blade_bnd2,blade_bnd3,blade_bnd4,blade_bnd5,blade_bnd6,blade_bnd7,blade_bnd8];
  for i=1:16
    l2(i)=length(find(blade_bnd==blade_bnd(1,i)));
  end
  pos=find(l2==1);
  for i=1:8
    side_bnd(1,i)=blade_bnd(pos(1,i));
    I(i)=postint(fem,'0.25*0.003[m]*(1.25[kg/m^3]*normv_acpr^2-p_t_acpr^2/(1.25[kg/m^3]*117649[m^2/s^2]))*(x*ny_acpn-y*nx_acpn)', ...
    'unit','kg*m^2/s^2','recover','off','dl',side_bnd(1,i),'edim',1);
  end
  M(k+1)=sum(I);
end
save moment M -ascii