************************************************************************
      SUBROUTINE sim_indat(t)
************************************************************************
      IMPLICIT NONE
C 
C REVISONS:
C --------
C  02/23/05 c.neumeyer Final modifications for RWM
C  11/11/03 t.gibney   Change calculation of "t".
C  07/28/03 t.gibney   Change from double-precision to single-precision.
C  10/25/00 t gibney   use parameters for NBRANCHES, NCOILS, etc.
C  03/14/00 t gibney   set imeas along with das:  imeas computed
C                       by acq program in new scheme PSC_1_nn
C  06/11/99 t gibney   if newCycle, set das(*)=0, return t=tsop
C  12/02/98 c neumeyer added fixes during PTP-ECS-034
C  11/17/98 c neumeyer added modifications for collect data
C  11/06/98 c neumeyer added TF & CHI link requirement
C  10/21/98 c neumeyer various post FDR fixes
C  10/04/98 c neumeyer various fixes in preparation for FDR
C  09/16/98 b davis    bug fix in build call. Calls to dtime.
C  09/09/98 c neumeyer
C
C Parameters:
      real*4 t		! <r> current time for simulation
C Common block "include"
      include 'rtcdef.inc'
C Local variables
      integer n,nc,ndtstep
      real*4 ibranch
C
C   simulation sub mode /////////////////////////////////////////
C
			!Simulate power supply response to current
			! outputs.  Afterwards, increment 't'
			!=================================================
      if (newCycle) then		!note: is set .FALSE. in main loop
         t = tsop
         nshot = nshot+1
         do n=1,NCHAN1R
            das(n) = 0.0
         enddo
         nsopeop=1
         return				!---> return if newCycle
      endif

      call sim_ps_response

      DO 280, ndtstep = 1,ndtsim
         CALL lr (nckt,nsim,ntype,ncs,nsc,mc,rc,m,r,x,vc,ic
     &         ,icdot,dtsim)
C
C   apply polarity constraints to currents
C   except for SPA which is inherently bipolar
C
         DO 260, n = 1,NCOILS
            IF (nap(n).NE.1.AND.n.le.NCOILS_V0) THEN
               IF (nplink(n).EQ.1.AND.ic(n).LT.0.) ic(n) = 0
               IF (nplink(n).EQ.-1.AND.ic(n).GT.0.) ic(n) = 0.
            ENDIF
260      CONTINUE

         IF (ndtstep.eq.ndtin) THEN
C
C   update input data during this step for next command calculation
C
C
C   update circuit resistances
C
            DO 265 n=1,NCOILS
               rc(n) = r20c(n)*(1.+coeff*(tcoil(n)-20))+ecr(n)
265         CONTINUE
C
C   update branch currents
C
            DO 270, n = 1,NBRANCHES
C
C   establish circuit number nc
C
                nc = ncb(n)
C
C   skip locked out circuits
C
                IF (nlock(nc).EQ.1) GOTO 270
C
C   skip branch if multiplier is set to zero
C
                IF (kpb(n).EQ.0) GOTO 270
C
C   branch current = simulated circuit current distributed amongst
C   parallels in accordance with scaling and link status
C
                ibranch = ic(nc)/np(nc)*abs(kpb(n))
C
C   check current direction vs. branch polarity constraint
C   except for SPA which is inherently bipolar
C
      if(n.gt.NCOILS_V0) GO TO 267
                IF (ic(nc)*nplink(nc).GE.0..AND.kpb(n).LT.0) ibranch=0.
                IF (ic(nc)*nplink(nc).LT.0..AND.kpb(n).GT.0) ibranch=0.
267   CONTINUE
C
C   scale simulated measurements
C
                imeas(nmeas(n,1)) = ibranch
                imeas(nmeas(n,2)) = ibranch
                das(nmeas(n,1)) = ibranch/kmeas(n,1)
                das(nmeas(n,2)) = ibranch/kmeas(n,2)
270         CONTINUE
        END IF
C
C   input data update complete
C
280   CONTINUE
C
C   simulation iterations complete
C
!!!      if(nsopeop.eq.1)t=t+dt
!!!      t = t+dt
      t = (nint(t/dt) + 1) * dt			! 11-Nov-2003  TRG
      return
      end



************************************************************************
      SUBROUTINE sim_ps_response
************************************************************************
      IMPLICIT NONE
C
C   02/6/02 c. neumeyer corrected voltages when one twin locked out
C Parameters
      ! <none>
C Include common block
      include 'rtcdef.inc'
      include 'pcs_common.inc'
C Local variables
      integer n,npolar
      REAL*4 anew
      REAL*4 pn,qn,vconv
C
C   vsim is modified version of vc which takes
C   the simulated (lumped) voltage and current results for twin
C   circuits under common control and includes in each of the twin
C   circuit entries 1/2 of the total voltage and the common current
C   Note: In practice, in case operation is not exactly symmetric
C   about the midplane, the twin voltages would not be exactly
C   equal due to different mutually coupled voltages
C

C
C   simulation updates //////////////////////////////////////////////
C
C
C   initialize v & i output arrays per simulation results
C
      DO 25 n=1,NCOILS
         vsim(n)=vc(n)
25    continue
C
C   modify v & i output arrays based on twin circuit usage
C
      DO 50, n = 1,NCOILS_V0
        IF (nlock(n).EQ.1) THEN
C
C   circuit is locked out
C
C   if twin circuit is active, divide voltage between twins
C   and set load currents equal
C
C  PF1aU and PF1aL
C
        IF (n.EQ.2.AND.nlock(3).NE.1) THEN
          vsim(2) = vc(3)/2.
          vsim(3) = vsim(2)
        ENDIF
        IF (n.EQ.3.AND.nlock(2).NE.1) THEN
          vsim(2) = vc(2)/2.
          vsim(3) = vsim(2)
        ENDIF
C
C  PF2U and PF2L
C
        IF (n.EQ.5.AND.nlock(6).NE.1) THEN
          vsim(5) = vc(6)/2.
          vsim(6) = vsim(5)
        ENDIF
        IF (n.EQ.6.AND.nlock(5).NE.1) THEN
          vsim(5) = vc(5)/2.
          vsim(6)=vsim(5)
        ENDIF
C
C  PF3U and PF3L
C
        IF (n.EQ.7.AND.nlock(8).NE.1) THEN
          vsim(7)=vc(8)/2.
          vsim(8)=vsim(7)
        ENDIF
        IF (n.EQ.8.AND.nlock(7).NE.1) THEN
          vsim(7) = vc(7)/2.
          vsim(8)=vsim(7)
        ENDIF
      END IF
50    CONTINUE

      DO 400, n = 1,NCOILS
         IF (nlock(n).EQ.1) GOTO 400
C
C   determine power supply voltage based on commands
C
C   voltage is zero if pos and neg converters are suppressed
C
         IF (ncvt1(n).EQ.0.AND.ncvt2(n).EQ.0) THEN
            vc(n) = 0.
            ntype(n) = 3
            GOTO 400
         ENDIF
C
C   determine which converter branch is active
C
         IF (nap(n).EQ.0) THEN
C
C   antiparallel not enabled, normal operation
C
C   constrain current to one direction
C
            IF (iload(n)*nplink(n).LT.0.) iload(n) = 0.
         ENDIF
C
C   limit rate of change of increase of effective alpha
C
         IF (n .LE. NCOILS_V0) then	!original Transrex pwr supplies
            IF (a1(n)-a1old(n).GT.damax) a1(n)= a1old(n)+damax
            IF (a2(n)-a2old(n).GT.damax) a2(n)= a2old(n)+damax
         ENDIF
         a1old(n) = a1(n)
         a2old(n) = a2(n)
C
C   determine which branch is active
C
         IF (iload(n)*nplink(n).GE.0.) THEN
C
C   positive converter active
C
            npolar = 1
            anew=a1(n)
         ELSE
C
C   negative converter active
C
            npolar = -1
            anew=a2(n)
         ENDIF

         if (n .LE. NCOILS_V0) then	!original Transrex coil set only
            if (n .ne. IDX_CHI .or. nctFlag) then   ! added to implement new CHI mode REH - 6/2k
               vconv = vmax(n)*cos(anew)-rps(n)*abs(iload(n))
            else
               if (anew .le. pi/3) then
                  vconv =  vmax(n)*cos(anew)-rps(n)*abs(iload(n))
               elseif (anew .gt. pi/3 .and. anew .le. 2.*pi/3) then
                  vconv = vmax(n)*(1.+sin(pi/6-anew))-
     $                    rps(n)*abs(iload(n))
               else
                  vconv = 0.
               endif
               if (vconv .lt. 0.) vconv = 0.
            endif
                
C
C   if inverting, limit invert voltage
C
            IF (vconv.LT.finv*vmax(n)) vconv = finv*vmax(n)
C
C   adjust polarity according to active converter and link status
C

         ELSE		!else SPA/RWM power supply ...
		anew=a1(n)
		npolar=1
            vconv = vmax(n) * anew/SPA_MAX_VREF
         ENDIF
         vc(n) = vconv*npolar*nplink(n)
400   CONTINUE
C
C   set voltage of twin circuits, compute composite power
C
      p = 0.
      q = 0.
      DO 450, n = 1,NCOILS
        IF (nlock(n).EQ.0) GOTO 440
C
C   circuit is locked out
C
C   if twin circuit is active, apply voltage to twin
C
C  PF1aU and PF1aL
C
        IF (n.EQ.2.AND.nlock(3).NE.1) THEN
          vc(3) = vc(3)*2.
          vc(2) = 0.
        ENDIF
        IF (n.EQ.3.AND.nlock(2).NE.1) THEN
          vc(2) = vc(2)*2.
          vc(3) = 0.
        ENDIF
C
C  PF2U and PF2L
C
        IF (n.EQ.5.AND.nlock(6).NE.1) THEN
          vc(6) = vc(6)*2.
          vc(5) = 0.
        ENDIF
        IF (n.EQ.6.AND.nlock(5).NE.1) THEN
          vc(5) = vc(5)*2.
          vc(6) = 0.
        ENDIF
C
C  PF3U and PF3L
C
        IF (n.EQ.7.AND.nlock(8).NE.1) THEN
          vc(8) = vc(8)*2.
          vc(7) = 0.
        ENDIF
        IF (n.EQ.8.AND.nlock(7).NE.1) THEN
          vc(7) = vc(7)*2.
          vc(8) = 0.
        ENDIF
440     CONTINUE
C
C   composite power
C
      IF (nfault(n,1).EQ.1) GOTO 450
        pn = vc(n)*ic(n)
        qn = sqrt((vmax(n)*ic(n))**2-pn**2)
        p = p+pn
        q = q+qn
450   CONTINUE
      s = sqrt(p**2+q**2)

      return
      end



C   ************************************************************
        SUBROUTINE lr
     &      (nckt,nsim,ntype,ncs,nsc,mc,rc,m,r,x,vc,ic,icdot,dt)
C   ************************************************************
C
C   simulate l-r circuit behavior
C
C   nckt is total number of circuits
C
C   nsim is number of active circuits
C
C   ntype gives ckt types
C
C   if ntype = 0, specification is v
C   if ntype = 1, specification is i
C   if ntype = 2, circuit is open
C   if ntype = 3, circuit is shorted
C   if ntype = 4, voltage feedback control
C   if ntype = 5, current feedback control
C
C   ncs contains total ckt number as function of sim circuit number
C   nsc contains sim circuit number as function of total circuit number
C
C   mc and rc are ind and res matrices of order nckt x nckt
C   l and r are ind and res matrices of order nsim x nsim
C
C
C   vc, ic, and icdot are of order nc
C   v, i, and idot are of order nsim
C
        IMPLICIT        REAL*4(a-m,o-z)
	include 'rtcconfig.inc'
        DIMENSION ntype(NCOILS),ncs(NCOILS),nsc(NCOILS)
        DIMENSION m(NCOILS,NCOILS),r(NCOILS),mc(NCOILS,NCOILS),
     &    rc(NCOILS,NCOILS),x(NCOILS,NCOILS),y(NCOILS)
        DIMENSION v(NCOILS),i(NCOILS),idot(NCOILS),vc(NCOILS),
     &    ic(NCOILS),icdot(NCOILS)
C
C   load circuit values into simulation slots
C
        DO 100, n1 = 1,nsim
          v(n1) = vc(ncs(n1))
          i(n1) = ic(ncs(n1))
          idot(n1) = icdot(ncs(n1))
100     CONTINUE
        DO 300, n1 = 1,nsim
          IF (ntype(ncs(n1)).EQ.1) THEN
C
C   current controlled circuit, calculate voltage
C
            v(n1) = r(n1)*i(n1)
            DO 200, n2 = 1,nsim
              v(n1) = v(n1)+m(n1,n2)*idot(n2)
200         CONTINUE
          ENDIF 
C
C   update v - r * i vector
C
          y(n1) = v(n1)-r(n1)*i(n1)
300     CONTINUE
        DO 500, n1 = 1,nsim
          idot(n1) = 0.
          DO 400, n2 = 1,nsim
            idot(n1) = idot(n1)+x(n1,n2)*y(n2)
400       CONTINUE
          i(n1) = i(n1)+idot(n1)*dt
500     CONTINUE 
C
C   finalize data
C
        DO 700, n1 = 1,NCOILS
          IF (ntype(n1).EQ.2) THEN
C
C   calculate induced voltages in open circuits
C
            vc(n1) = 0.
            DO 600, n2 = 1,NCOILS
              IF (ntype(n2).NE.2) THEN
                vc(n1) = vc(n1)+mc(n1,n2)*idot(nsc(n2))
              ENDIF
600         CONTINUE
            ic(n1) = 0.
            icdot(n1) = 0.
          ELSE 
C
C   transfer sim results to total circuit matrices
C
            vc(n1) = v(nsc(n1))
C     IF(vc(n1).gt.10.e3.or.vc(n1).lt.-10.e3)THEN
C   write(6,*)'voltage limiting ckt.',n1
C   IF(vc(n1).gt.0.)vc(n1)=10.e3
C   IF(vc(n1).lt.0.)vc(n1)=-10.e3
C     END IF
            ic(n1) = i(nsc(n1))
            IF (ic(n1).GT.2.e6.OR.ic(n1).LT.-2.e6) THEN
              WRITE(6,*)'current beyond limit ckt.',n1
              STOP
            ENDIF
            icdot(n1) = idot(nsc(n1))
          ENDIF
700     CONTINUE
        RETURN
        END 
