/out,m_induc,out
/com,m_induc.mac
!  Computes the L matrix vs current for 2D models
!
!   arg1 = winding file name (no default)
!   arg2 = electrical angle for the current specification (defauts to 90)
!   arg3 = position of the rotor (mechanical angle measured
!          in degrees CCW from +X axis)  (defaults to zero)
!   arg4 = starting current (not amp-turns) (A) (defaults to 1 amp-turn)
!   arg5 = ending current (not amp-turns)   (A)
!          defaults to the starting current for a single solution
!   arg6 = number of total solutions (defaults to one)
!
!
alls
*get,ecnt,elem,,count
*if,ecnt,le,0,then
  /out
  /nerr
  *msg,error
  There are no elements in the model.
*endif
*get,_mcpu1,active,,time,cpu
/nerr,0,1e5
*if,ecnt,le,0,:end
h_mh=1000   !   converts H to mH
_arg1=

*if,arg9,eq,1,then
  *ask,_arg1,Winding file name
  *ask,_arg2,Electrical angle for current,90
  *ask,_arg3,Initial position of rotor,0
  *ask,_arg4,Initial current,1
  *ask,_arg5,Ending current,1
  *ask,_arg6,Number of solutions,1
  *ask,go_now,1 to continue OR 0 to Stop,0
  arg1=_arg1
  arg2=_arg2
  arg3=_arg3
  arg4=_arg4
  arg5=_arg5
  arg6=_arg6
*else
 go_now=1
*endif
*msg,info,arg1,arg2,arg3,arg4,arg5,arg6
arg1,arg2,arg3,arg4,arg5,arg6:  %c  %g %g %g %g %g %g

*if,arg6,lt,0,:end

*if,go_now,eq,0,:end

esel,none
cmsel,,stator
*get,_elems,elem,,count

esel,none
cmsel,,rotor
*get,_elemr,elem,,count
_err=0
*if,_elems*_elemr,le,0,then
  _err=1
  /nerr
  /out
  *msg,error  
  Part of the model is missing-no action
*endif
/nerr,0,1e5
*if,_err,gt,0,:end



/com, build the winding data
_wfnam=arg1
*if,_nomsg,eq,1,then
  _mg1=1
*endif
*stat,_wfnam
b2wndsc,_wfnam,'s_coil'
/out,m_induc,out,,append
*if,_err,eq,1,:end

*if,_nomsg,eq,1,then
  _mg1=1
*endif
/com, move / connect the rotor
mvrotor,arg3,1
/out,m_induc,out,,append

!*if,arg2,le,0,then
!  arg2=90
!*endif

*if,arg4,le,0,then
  arg4=1
*endif

*if,arg5,le,0,then
 arg5=arg4
*endif

*if,arg6,le,0,then
  arg6=1
*endif

_nst=arg6

*msg,info,arg1,arg2,arg3,arg4,arg5,arg6
arg1,arg2,arg3,arg4,arg5,arg6:  %c  %g %g %g %g %g %g

*set,m_induc
*dim,m_induc,table,_nst,6
m_induc(0,0)=1e-6

*if,_nst,gt,1,then
 *do,_i1,1,_nst
   m_induc(0,_i1)=_i1
 *enddo
*endif


*if,arg4,gt,arg5,then
 _t=arg5
 arg5=arg4
 arg4=_t
 _t=
*endif

*if,_nst,eq,1,then
 _delcur=0
*else
 _delcur=(arg5-arg4)/(_nst-1)
*endif
*if,_delcur,lt,.001,then
  _nst=1
*endif

*stat,_delcur

*if,arg2,eq,0,then
  arg2=.001
*endif
*if,arg2,eq,180,then
  arg2=179.999
*endif
*if,stkthk,le,0,then
  _stk=1
*else
  _stk=stkthk
*endif
_nindt=(_mx_ph+1)*_mx_ph/2
_ttmax=-1
*do,_im1,1,_nst
 /gopr
 *msg,info,_im1
 Before the call: cycle _im1: %i
 _ccurr=arg4+_delcur*(_im1-1)
 _mg1=1
 *stat,_ccurr
 *msg,info,arg2
 arg2:  %g

 s_ind,arg2,_ccurr
 /out,m_induc,out,,append
 _fac=ggeom*_stk*h_mh
 *msg,info, _fac,ggeom,_stk,h_mh
  _fac,ggeom,_stk,h_mh:  %g  %g  %g  %g


 m_induc(_im1,1)=induc(1,1)*_fac
 m_induc(_im1,2)=induc(2,2)*_fac
 m_induc(_im1,3)=induc(3,3)*_fac
 m_induc(_im1,4)=induc(1,2)*_fac
 m_induc(_im1,5)=induc(1,3)*_fac
 m_induc(_im1,6)=induc(2,3)*_fac
 m_induc(_im1,0)=_ccurr
 *vabs,,1
 *vscfun,_tmax,max,induc(1,1)
 *vabs
 _tmax=_tmax*_fac
 *if,_tmax,gt,_ttmax,then
   _ttmax=_tmax
 *endif
*enddo
*get,_mcpu2,active,,time,cpu
mdcpu=_mcpu2-_mcpu1
*abbr,_vwrt1,*vwrit,m_induc(1,0),m_induc(1,1),m_induc(1,2)
/nopr
/out,m_induc,sum
*msg,info
______________MATRIX SUMMARY CALCULATION____________________
*msg,info,arg1
 Winding file:___________ %c
*msg,info,arg2
 Electrical angle:_______ %g
*msg,info,arg3
 Rotor position:_________ %g (CCW from +X axis)
*msg,info,facsym
Symmetry factor:_________ %i
*msg,info,stkthk
Stack length (inches):___ %g
*msg,info

*msg,info
__Current_______________Inductance (mH) Matrix_______________
*msg,info
__(A)_______Laa______Lbb______Lcc______Lab______Lac______Lbc_
!234567890!234567890!234567890!234567890!234567890!234567890
!11111222222223333333
*msg,info

*if,_ttmax,lt,100,then
_vwrt1,m_induc(1,3),m_induc(1,4),m_induc(1,5),m_induc(1,6)
(f6.2,3x,6(1x,f8.3))
*else
_vwrt1,m_induc(1,3),m_induc(1,4),m_induc(1,5),m_induc(1,6)
(f6.2,3x,6(1x,f8.1))
*endif
*msg,info
____________________________________________________________
*vwrit,mdcpu
('Total CPU time (Sec): ',f6.0)
/out
*if,_nomsg,eq,0,then
 *uili,m_induc,sum
*endif

*abbr,_vwrt1
:end
_mg1=0
/out