/out,pmstat,out
/com, pm_stat.mac    11/05/01   added air/housing outside stator
/com,pm_stat.mac     12/25/00   original version
!   *abbr,ames,stat
!   pmstat.mac
!   Generates the PM stator based on parameters
!    arg1 = mesh refinement, 1 is minimum number of elements, 5 is most 
!           most refinement
!    arg2 = 0  use lower order element for the iron
!         > 0  use higher order element for the iron
!    arg3 = number of stator teeth to be generated
!         = 0 (default), no action, no additional teeth are generated
!         ne 0 Build arg3 number of teeth
!    arg4 = 0  use lower order element for the magnet and air gap
!         > 0  use higher order element for the magnet and air gap
!
!

fini
/prep7
/uis,msgpop,3
!  /nerr,0,-1

_error=0
/nerr,1,1e6
*get,_mod_id,parm,stat_id,type
*if,_mod_id,eq,3,then
*if,stat_id,ne,'pmstat',then
 /nerr,1,1e6
 /out
 *msg,error,stat_id
 The model type ( %c ) does not correspond to a Permanent magnet stator
 /out,pmstat,out,,append
 _error=1
*endif
*else
 /out
 _error=1
*msg,error
 The model type identifier specified is not a "character" type parameter&
 or it was not specified
 /out,pmstat,out,,append
*endif


*if,_error,eq,1,:err

/nerr,1,-1
!  geometry error check

_error=0

_nm_vn=4
*set,_nm_vv
*dim,_nm_vv,,_nm_vn
_nm_vv(1)=r4,r5,nsp,gap
*set,_nm_v
*dim,_nm_v,char,_nm_vn
_nm_v(1)='r4','r5','nsp','gap'
*do,_i,1,_nm_vn
*if,_nm_vv(_i),eq,0,then
 *if,_i,ne,6,then
  *msg,error,_nm_v(_i)
  Parameter { %c }  is negative or zero.
  _error=_error+1
 *endif
*endif
*enddo


*if,src2+src1,le,0,then
  *if,r3,le,0,then
   *msg,error
   Parameter {rr2} OR {magh} are negative or zero.
  _error=_error+1
  *endif
*endif

*if,smagl+smagw,eq,0,then
  *msg,error
  Parameter {smagl} and {smagw} are negative or zero.
  _error=_error+1
*endif
*if,smagw,gt,180,then
  *msg,error
  The magnet angle width cannot be greater than 180
  _error=_error+1
*endif
/nerr,0,1e6

*if,_error,eq,1,:err
/nerr,0,1e4
*if,arg1,eq,0,then
 arg1=2
*endif
_cd1=1    !   radial divisons for the air space 
!             between two phases in the same slot
*get,_elmtyp,etyp,,num,max
*if,_elmtyp,eq,0,then
 _elmtyp=1
*endif
_elctyp=_elmtyp+1

*if,arg2,eq,0,then
  et,_elmtyp,13
*else
  et,_elmtyp,53
*endif

*if,arg4,eq,0,then
  et,_elctyp,13
*else
  et,_elctyp,53
*endif

*if,arg6,eq,1,then
  et,_elctyp,53,3
*endif



immed
esel,none
nsel,none
asel,none
ksel,none
lsel,none
cmsel,,stator
cpdel,all,,,any
cedel,all,,,any
acle,all
adel,all,,,1
ldel,all,,,1
kdel,all
edel,all
ndel,all
! numc,kp
! numc,line
! numc,area
numc,elem
numc,node
numc,cp
numc,ce

alls
cm,_ckp,kp
cm,_cls,line
cm,_car,area
cmgrp,_cmp,_ckp,_cls,_car




!   determine the current maximum real constant set
esel,all
*get,_elmx,elem,,num,max
_rl_strt=0
*if,_elmx,ne,0,then
*set,_mskv
*dim,_mskv,,_elmx
*vget,_mskv(1),elem,1,esel
*set,_r_e
*dim,_r_e,,_elmx
*vmask,_mskv(1)
*vget,_r_e(1),elem,1,attr,real
*vmask,_mskv(1)
*vscfun,_rl_strt,max,_r_e(1)
_zz=100*nint(_rl_strt/100)
*if,_zz,lt,_rl_strt,then
  _rl_strt=_zz+51
*else
  _rl_strt=_zz
*endif
_zz=
*set,_mskv
*set,_r_e
esel,none
*endif



!  compute locations
_mxkp=11          !  11/06/01
*if,t_air,gt,0,then
  _mxkp=13
*endif
*if,t_hous,gt,0,then
  _mxkp=_mxkp+12
*endif

*set,_kp_n_
*dim,_kp_n,,_mxkp
_60=180
*set,_x
*dim,_x,,_mxkp
*set,_y
*dim,_y,,_mxkp
*afun,deg

*if,rbc,le,0,then
 *if,src1+src2,eq,0,then
  rbc=r3-gap/2
 *else
  rbc=src1+src2-gap/2
 *endif
*endif

*if,nsp,eq,0,then
 _err=2
 nsp=1
*endif

_x(1)=rbc
_y(1)=0
_x(2)=rbc*cos(180/nsp)
_y(2)=rbc*sin(180/nsp)

*if,src1+src2,eq,0,then
 _x(3)=r3
 _y(3)=0
*else
 _x(3)=src1+src2
 _y(3)=0
*endif



*if,src1+src2,ne,0,then
 ! the shaped pole option
 *if,smagw,ne,0,then
  !  magnet width is by the angle 
  t_b=tan(smagw/nsp)
  _x(4)=src1+sqrt(src1**2-(1+t_b**2)*(src1**2-src2**2))
  _x(4)=_x(4)/(1+t_b**2)
  _y(4)=_x(4)*t_b
 *elseif,smagl,ne,0,then
  _x(4)=src1+sqrt(src2**2-(smagl/2)**2)
  _y(4)=smagl/2
 *else
  _err=3
  _x(4)=0
  _y(4)=0
 *endif
!  *elseif,smagh,ne,0,then
*elseif,r3,ne,0,then
 ! the constant radius option
 *if,smagw,ne,0,then
  !  magnet width is by the angle 
  _x(4)=(r3)*cos(smagw/nsp)
  _y(4)=(r3)*sin(smagw/nsp)
 *elseif,smagl,ne,0,then
  _x(4)=sqrt((r3)**2-(smagl/2)**2)
  _y(4)=smagl/2
 *else
  _err=3
  _x(4)=0
  _y(4)=0
 *endif
*else
  _err=3
  _x(4)=0
  _y(4)=0
*endif

_x(5)=r4
_y(5)=0

*if,sloaf,eq,1,then
 _x(6)=sqrt(r4**2-_y(4)**2)
 _y(6)=_y(4)
*else
 _x(6)=r4*cos(smagw/nsp)
 _y(6)=r4*sin(smagw/nsp)
*endif
!  
_x(7)=(r4)*cos(_60/nsp)
_y(7)=(r4)*sin(_60/nsp)

_x(8)=r5
_y(8)=0

_x(9)=(r5)*cos(_60/nsp)
_y(9)=(r5)*sin(_60/nsp)

*if,src1+src2,ne,0,then
 _x(10)=src1
 _y(10)=0
*else
 _x(10)=0
 _y(10)=0
*endif

_r_1=sqrt(_x(4)**2+_y(4)**2)
_x(11)=_r_1*cos(_60/nsp)
_y(11)=_r_1*sin(_60/nsp)

*if,t_air,gt,0,then
 _x(12)=r5+t_air
 _y(12)=0
 _x(13)=(r5+t_air)*cos(_60/nsp)
 _y(13)=(r5+t_air)*sin(_60/nsp)
*else
 t_air=0
*endif

*if,t_hous,gt,0,then
 _x(14)=r5+t_air+t_hous
 _y(14)=0
 _x(15)=_x(14)*cos(_60/nsp)
 _y(15)=_x(14)*sin(_60/nsp)
*endif

csys
ksel,none
*do,_ikp,1,11  ! _mxkp    !  11/05/01
/gopr
 k,,_x(_ikp),_y(_ikp)
 _kp_n(_ikp)=kp(_x(_ikp),_y(_ikp),0)
*enddo

*if,t_air,gt,0,then         !  11/05/01   start
 *do,_ikp,12,13
  k,,_x(_ikp),_y(_ikp)
  _kp_n(_ikp)=kp(_x(_ikp),_y(_ikp),0)
 *enddo
*endif

*if,t_hous,gt,0,then
 *do,_ikp,14,15
  k,,_x(_ikp),_y(_ikp)
  _kp_n(_ikp)=kp(_x(_ikp),_y(_ikp),0)
 *enddo
*endif                     !  11/05/01   end

cm,sta_k,kp

!  the  air gap next to the magnet
asel,none
lsel,none
csys,1
l,_kp_n(1),_kp_n(2)
l,_kp_n(2),_kp_n(11)
l,_kp_n(11),_kp_n(4)
*if,src1+src2,ne,0,then
 larc,_kp_n(4),_kp_n(3),_kp_n(10),src2
*else
 csys,1
 l,_kp_n(4),_kp_n(3)
 csys
*endif
l,_kp_n(3),_kp_n(1)
a,_kp_n(1),_kp_n(2),_kp_n(11),_kp_n(4),_kp_n(3),
*get,_ar_1,area,,num,max
aatt,1,6,_elctyp
ksel,,,,_kp_n(3)
ksel,a,,,_kp_n(4)
ksel,a,,,_kp_n(11)
lslk,,1
*get,_lsmx,line,,num,max
*get,_lsmn,line,,num,min
lccat,_lsmn,_lsmx
cm,sta_a,area
!   the air beside the magnet
ksel,all
asel,none
lsel,none
csys,1
l,_kp_n(11),_kp_n(7)
l,_kp_n(7),_kp_n(6)
csys
l,_kp_n(6),_kp_n(4)
a,_kp_n(4),_kp_n(11),_kp_n(7),_kp_n(6),
aatt,1,6,_elmtyp
cmsel,a,sta_a
cm,sta_a,area

!  the magnet 
! locate the next available CS, above 50
! *do,_i1,50,100                ! removed 2/15/96
!  *get,_esys,cdsy,_i1,attr,kcs ! removed 2/15/96
!  *if,_esys,eq,-1,exit         ! removed 2/15/96
! *enddo                        ! removed 2/15/96
! _mag_esy=_i1                  ! removed 2/15/96
! _esys=     $   _i1=           ! removed 2/15/96

*get,_mag_esy,cdsy,,num,max     !  added 2/15/96
*if,_mag_esy,lt,50,then         !  added 2/15/96
  _mag_esy=50                   !  added 2/15/96
*else                           !  added 2/15/96
  _mag_esy=_mag_esy+1           !  added 2/15/96
*endif                          !  added 2/15/96
!   local,_mag_esy,0            !  removed 6/26/96
local,_mag_esy,rad_mag          !  6/26/96
csys                            !  6/26/96


 asel,none
 lsel,none
 ksel,,,,_kp_n(6)
 ksel,a,,,_kp_n(4)
 ksel,a,,,_kp_n(3)
 lslk,,1
 ksel,all
 l,_kp_n(3),_kp_n(5)
 csys,1
 l,_kp_n(5),_kp_n(6)
 al,all
 aatt,6,3,_elctyp,_mag_esy
 cmsel,a,sta_a
 cm,sta_a,area

 !  the iron return path for the rotor
 asel,none
 lsel,none
 ksel,,,,_kp_n(5)
 ksel,a,,,_kp_n(6)
 ksel,a,,,_kp_n(7)
 lslk,,1
 ksel,all
 csys,1
 l,_kp_n(5),_kp_n(8)
 l,_kp_n(8),_kp_n(9)
 l,_kp_n(9),_kp_n(7)
 al,all

aatt,5,2,_elmtyp
cmsel,a,sta_a
cm,sta_a,area


/com, air gap outside the stator   !  11/05/01   start
*if,t_air,gt,0,then
  asel,none
  a,_kp_n(8),_kp_n(12),_kp_n(13),_kp_n(9)
  aatt,1,10,1
  cmsel,a,sta_a
  cm,sta_a,area
*endif


*if,t_hous,gt,0,then
  asel,none
  *if,t_air,eq,0,then
   a,_kp_n(8),_kp_n(14),_kp_n(15),_kp_n(9)
  *else
   a,_kp_n(12),_kp_n(14),_kp_n(15),_kp_n(13)
  *endif
  aatt,7,10,1
  cmsel,a,sta_a
  cm,sta_a,area
*endif                         !  11/05/01  end



!   set meshing divisions
!   number divisions in the radial gap
*set,_gapd
*dim,_gapd,,5
_gapd(1)=2,2,3,4,5
!   aspect ratio of tooth along the tooth face
*set,_asprg
*dim,_asprg,,5
_asprg(1)=2,2,1,1,1
!  factor applied to kesize for keypoints in the air gap
!  based on the gap/_gapd(arg1) value
*set,_kpz1
*dim,_kpz1,,5
_kpz1(1)=4,3,2,1,1
!  factor applied to kesize for keypoints at the back of the magnet 
*set,_kpz2
*dim,_kpz2,,5
_kpz2(1)=24,18,9,6,3
!  factor applied to kesize for keypoints at the back of the iron
*set,_kpz3
*dim,_kpz3,,5
_kpz3(1)=24,20,18,6,4

ksel,,,,_kp_n(1)
ksel,a,,,_kp_n(3)
lslk,,1
lesi,all,,,_gapd(arg1)
ksel,,,,_kp_n(11)
ksel,a,,,_kp_n(2)
lslk,,1
lesi,all,,,_gapd(arg1)

ksel,,,,_kp_n(1)
ksel,a,,,_kp_n(3)
ksel,a,,,_kp_n(2)
ksel,a,,,_kp_n(4)
ksel,a,,,_kp_n(11)

lesi,all,_asprg(arg1)*gap/_gapd(arg1)/2

*if,t_air,gt,0,then       !  11/05/01  start
 ksel,,,,_kp_n(8)
 ksel,,,,_kp_n(12)
 lslk,,1
 lesi,all,,,_gapd(arg1)
 ksel,,,,_kp_n(7)
 ksel,,,,_kp_n(13)
 lslk,,1
 lesi,all,,,_gapd(arg1)
*endif


*if,t_hous,gt,0,then
 ksel,,,,_kp_n(14)
 ksel,,,,_kp_n(12)
 lslk,,1
 lesi,all,,,_gapd(arg1)*2
 ksel,,,,_kp_n(15)
 ksel,,,,_kp_n(13)
 lslk,,1
 lesi,all,,,_gapd(arg1)*2
*endif                         !  11/05/01  end


_gap_d=gap/_gapd(arg1)/2

ksel,all
kesi,_kp_n(3),_kpz1(arg1)*_gap_d
kesi,_kp_n(4),_kpz1(arg1)*_gap_d
kesi,_kp_n(11),_kpz1(arg1)*_gap_d

kesi,_kp_n(5),_kpz2(arg1)*_gap_d
kesi,_kp_n(6),_kpz2(arg1)*_gap_d
kesi,_kp_n(7),_kpz2(arg1)*_gap_d

kesi,_kp_n(8),_kpz3(arg1)*_gap_d
kesi,_kp_n(9),_kpz3(arg1)*_gap_d

cmsel,,sta_a
esha
!   mesh the lower order elements first
cmsel,,sta_a
*if,arg2,eq,0,then
  asel,r,type,,_elmtyp
  ames,all
*endif
cmsel,,sta_a
*if,arg4,eq,0,then
  esha,2
  ames,_ar_1
  esha
  asel,r,type,,_elctyp
  ames,all
*else
  esha,2
  ames,_ar_1
  esha
*endif
cmsel,,sta_a
ames,all

cmsel,,sta_a
esla
nsle
lsla
ksll
cm,sta_a,area
cm,sta_l,line
cm,sta_k,kp
cm,sta_e,elem
cm,sta_n,node
cmgrp,stator,sta_a,sta_l,sta_k,sta_e,sta_n
dsys
esel,r,mat,,6
cmsel,,stator

*if,arg3,ne,0,then
!          1   2 3 4 5    6     7   8
  polegen,arg3,0,0,0,0,'stator',2,ggeom
  /out,pmstat,out,,append
  csys,0
  mrk_nod=node(_x(3)*ggeom,_y(3)*ggeom,0)    !  marker node used in the 
                                             !  rotor movement
*endif
!   apply az at outer nodes
cmsel,,stator
csys,1
nsel,,ext
*get,_rdmx,node,,mxloc,x
nsel,r,loc,x,_rdmx-.0001,_rdmx+.0001
d,all,az
/pbc,a,1
cmsel,,stator



cmsel,,stator

*if,rad_mag,eq,1,then
*get,_mzz,mgxx,6,temp,20
mp,mgxx,7,_mzz*(-1)
*get,_mzz,murx,6,temp,20
mp,murx,7,_mzz
_zz=
!   adjust the magnets to alternate in-out
_iflip=1
 *do,_i1,1,nsp
 /gopr
 *if,_iflip,eq,-1,then
   cmsel,,stator
   esel,r,mat,,6
   nsle
   csys,1
   nsel,r,loc,y,360/nsp*((2*_i1-1)/2-1)-.0001,360/nsp*((2*_i1-1)/2)+.0001
   esln,,1
   emod,all,mat,7
   _iflip=_iflip*(-1)
 *else
   _iflip=_iflip*(-1)
 *endif
 *enddo
*endif



/pnum,mat,1
csys
*get,_xmn,node,,mnloc,x
*get,_xmx,node,,mxloc,x
*get,_ymn,node,,mnloc,y
*get,_ymx,node,,mxloc,y
_dify=_ymx-_ymn
_difx=_xmx-_xmn
_eify=_ymx+_ymn
_eifx=_xmx+_xmn
_rmx=_dify
*if,_difx,gt,_dify,then
   _rmx=_difx
*endif
 /foc,1,_eifx/2,_eify/2
 /dis,1,1.05*_rmx/2
!  clean out the solid model
alls
cmsel,u,_cmp
! adel,all,,,1
!ldel,all,,,1
!kdel,all
!numc,kp
!numc,line
!numc,area
/auto
/num,1
eplo
immed,1


:err

/out