      implicit real*8 (a-m,o-z)
c
c   nc is max number of circuits
c
      parameter(nc=75)
c
c   np is max number of breakpoints for each circuit
c
c   no is max number of output v & i
c
      parameter(np=200,no=20)
      parameter(nump=5)
c
c   nump is max no. of points at which Bz vs. time is to be computed
c
      dimension lc(nc,nc),rc(nc),l(nc,nc),r(nc),x(nc,nc)
      dimension ntype(nc),ncs(nc),nsc(nc),nsum(nc)
      dimension vc(nc),ic(nc),icdot(nc)
      dimension csac(nc)
      dimension ndat(nc)
      dimension d(nc,np,2),nd(nc)
      dimension noc(no),isum(no)
      dimension ncomp(10),wcomp(10),icomp(10),dbdi(nc,nump),bz(nump)
      dimension nfc(nc),gp(nc),gi(nc),gd(nc),vlim(nc,2)
      dimension rps(nc),ecr(nc),error(nc),ierror(nc)
      dimension bpoint(np,2)
      character*12 elenam(nc)
      open(10,file='sim.d',status='old')
      open(11,file='nam.d',status='old')
      open(20,file='sim.v',status='old')
      open(30,file='sim.i',status='old')
      open(40,file='sim.b',status='old')
      open(50,file='bp.o',status='old')
      open(60,file='soi.o',status='old')
      open(70,file='eoi.o',status='old')
      open(80,file='csa.d',status='old')
      open(90,file='sim.j',status='old')
      open(95,file='sim.s',status='old')
c
c   read basic simulation data
c
c   nckt = total number of circuits
c   dt = time step
c   tstop = stop time
c   ndout = number of output data points
c
      read(10,*)nckt
      read(10,*)dt
      read(10,*)tstop
      read(10,*)ndout
      nskip=tstop/dt/ndout
      write(6,*)'sim.d timing data read completed'
c
c   read circuit number, circuit type, number of breakpoints
c     read circuit names
c
      d=0.
      do 200 n=1,nckt
      read(10,*)ncn,ntype(n),ndat(n),nsum(n)
      read(11,*)elenam(n)
c
c   read time/data breakpoints
c
      do 100 n1=1,ndat(n)
      read(10,*)d(n,n1,1),d(n,n1,2)
100   continue
200   continue
      write(6,*)'sim.d circuit data read completed'
c
c   nockt = number of output data circuits
c
      read(10,*)nockt,nsckt
c
c   read circuit numbers of output data circuits
c
      do 250 n=1,nockt
      read(10,*)noc(n)
250   continue
      write(20,888)(elenam(noc(n1)),n1=1,nockt)
      write(30,888)(elenam(noc(n1)),n1=1,nockt)
      write(6,*)'sim.d output spec read completed'
c
c   read feedback circuit numbers and settings
c
c   set voltage limit according to number of series pss
c   limit invert voltage
c   set rps to ps equivalent resistance
c
c   vdo and rpss per PPPL D-site power supplies
c
      vdo=1012.85
      kinv=-2./3.
      rpss=14.6e-3
      rps=0.
      read(10,*)nfbc
      if(nfbc.eq.0)go to 253
      do 252 n=1,nfbc
      read(10,*)nfc(n),gprop,gint,gder,nspss,nppss
      gp(nfc(n))=gprop
      gi(nfc(n))=gint
      gd(nfc(n))=gder
      vlim(nfc(n),1)=vdo*nspss
      vlim(nfc(n),2)=kinv*vdo*nspss
      rps(nfc(n))=rpss*nspss/nppss
252   continue
      write(6,*)'sim.d feedback spec read completed'
253   continue
c
c   read compensation circuit numbers and weighting
c
      read(10,*)ncc
      if(ncc.eq.0)go to 255
      do 254 n=1,ncc
      read(10,*)ncomp(n),wcomp(n)
254   continue
      write(6,*)'sim.d compensation spec read completed'
255   continue
c
c   if active compensation, read waveform for Bz at null target
c   and convert from gauss to tesla
c
		read(10,*)nbpoint
      if(nbpoint.eq.0)go to 257
		do 256 n=1,nbpoint
		read(10,*)bpoint(n,1),bpoint(n,2)
		bpoint(n,2)=bpoint(n,2)/1.e4
256		continue
      write(6,*)'sim.d Bz waveform read completed'
257   continue
      read(10,*)tsoi,teoi,tcomp
      write(6,*)'sim.d read completed'
c
c   read dbdi per amp data and cross sectional area data
c
      read(50,*)npoint
      do 260 n=1,nckt
      read(50,*)n1,(dbdi(n,n2),n2=1,npoint)
      read(80,*)n1,csac(n)
260   continue
      write(6,*)'b.o read completed'
      call build(nckt,ntype,nsim,lc,rc,l,r,x,ncs,nsc,ecr,rps)
      write(6,*)'build completed'
      nd=1
      nsoi=0
      neoi=0
      t=0.
      vc=0.
      ic=0.
      icdot=0.
      icomp=0.
      ncount=1
      ndb=1
c
c   enter time step loop
c
      do 900 n=1,tstop/dt
      if(ncc.ne.0)then
c
c   linear interpolation of Bz waveform for null target
c
      tlast=bpoint(ndb,1)
      tnext=bpoint(ndb+1,1)
      if(tnext.eq.tlast)then
		write(6,*)'Error: Tnext=Tlast In Bz Waveform'
		stop
      end if
c
c   if time has passed breakpoint, increment counter
c
      if(t.gt.tnext)then
		tlast=tnext
		ndb=ndb+1
		tnext=bpoint(ndb+1,1)
      end if
      dlast=bpoint(ndb,2)
      dnext=bpoint(ndb+1,2)
      bzref=dlast+(dnext-dlast)/(tnext-tlast)*(t-tlast)
      end if
      call fnull(nckt,ncc,ncomp,wcomp,icomp,dbdi,ic,bz,npoint,bzref)
      do 400 n1=1,nsim
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
      if(ntype(ncs(n1)).eq.2.or.ntype(ncs(n1)).eq.3)go to 400
c
c   henceforth ntype must be 0 or 1 or 4 or 5
c
c   linear interpolation
c
      tlast=d(ncs(n1),nd(ncs(n1)),1)
      tnext=d(ncs(n1),nd(ncs(n1))+1,1)
      if(tnext.eq.tlast)then
		write(6,*)'Error: Tnext=Tlast Ckt#',ncs(n1)
		stop
      end if
c
c   if time has passed breakpoint, increment counter
c
      if(t.gt.tnext)then
		tlast=tnext
		nd(ncs(n1))=nd(ncs(n1))+1
		tnext=d(ncs(n1),nd(ncs(n1))+1,1)
      end if
      dlast=d(ncs(n1),nd(ncs(n1)),2)
      dnext=d(ncs(n1),nd(ncs(n1))+1,2)
      data=dlast+(dnext-dlast)/(tnext-tlast)*(t-tlast)
      if(ntype(ncs(n1)).eq.0)then
c
c   voltage specification
c
      vc(ncs(n1))=data
      go to 400
      end if
      if(ntype(ncs(n1)).eq.1)then
c
c   current specification
c
      icdot(ncs(n1))=(dnext-dlast)/(tnext-tlast)
      go to 400
      end if
      if(ntype(ncs(n1)).eq.4)then
c
c   feedback controlled voltage
c
      tlimit=1./60./4.
      vlast=vc(ncs(n1))
      oerror=error(ncs(n1))
      error(ncs(n1))=data-vlast
      derror=(error(ncs(n1))-oerror)/dt
      ierror(ncs(n1))=ierror(ncs(n1))+error(ncs(n1))*dt
      vnext=gp(ncs(n1))*error(ncs(n1))+gi(ncs(n1))*ierror(ncs(n1))
     &+gd(ncs(n1))*derror
c
c   voltage rate limit
c
      dv=abs(vnext-vlast)
      dvlimit=vlim(ncs(n1),1)*dt/tlimit
      if(dv.gt.dvlimit)then
		dv=dvlimit
		if(vnext.gt.vlast)vnext=vlast+dv
		if(vnext.lt.vlast)vnext=vlast-dv
      end if
c
c   voltage clamp
c
      if(ic(ncs(n1)).ge.0.)then
		vmax=vlim(ncs(n1),1)
		vmin=vlim(ncs(n1),2)
	else
		vmax=-1.*vlim(ncs(n1),2)
		vmin=-1.*vlim(ncs(n1),1)
      end if
      if(vnext.gt.vmax)vnext=vmax
      if(vnext.lt.vmin)vnext=vmin
      vc(ncs(n1))=vnext
      go to 400
      end if
      if(ntype(ncs(n1)).eq.5)then
c
c   feedback controlled current
c
      tlimit=1./60./4.
c
c   set data to icomp if compensation circuit and after soi
c
      do 350 n2=1,ncc
      if(n1.eq.nsc(ncomp(n2)).and.t.ge.tcomp)data=icomp(n2)
350   continue
      vlast=vc(ncs(n1))
      oerror=error(ncs(n1))
      error(ncs(n1))=data-ic(ncs(n1))
      derror=(error(ncs(n1))-oerror)/dt
      ierror(ncs(n1))=ierror(ncs(n1))+error(ncs(n1))*dt
      vnext=gp(ncs(n1))*error(ncs(n1))+gi(ncs(n1))*ierror(ncs(n1))
     &+gd(ncs(n1))*derror
c
c   voltage rate limit
c
      dv=abs(vnext-vlast)
      dvlimit=vlim(ncs(n1),1)*dt/tlimit
      if(dv.gt.dvlimit)then
		dv=dvlimit
		if(vnext.gt.vlast)vnext=vlast+dv
		if(vnext.lt.vlast)vnext=vlast-dv
      end if
c
c   voltage clamp
c
      if(ic(ncs(n1)).ge.0.)then
		vmax=vlim(ncs(n1),1)
		vmin=vlim(ncs(n1),2)
	else
		vmax=-1.*vlim(ncs(n1),2)
		vmin=-1.*vlim(ncs(n1),1)
      end if
      if(vnext.gt.vmax)vnext=vmax
      if(vnext.lt.vmin)vnext=vmin
      vc(ncs(n1))=vnext
      go to 400
      end if
400   continue
      call lr(nckt,nsim,ntype,ncs,nsc,lc,rc,l,r,x,vc,ic,icdot,dt)
c
c   place complete current vector in output file if soi or eoi
c
997   format(i4,1e12.4)
      if(t.gt.tsoi.and.nsoi.eq.0)then
		do 990 n1=1,nckt
		write(60,997)n1,ic(n1)
990   continue
		nsoi=1
      end if
      if(t.gt.teoi.and.neoi.eq.0)then
		do 991 n1=1,nckt
		write(70,997)n1,ic(n1)
991   continue
		neoi=1
      end if
      t=t+dt
      ncount=ncount+1
      if(ncount.gt.nskip)then
      ncount=1
c
c   compute and store current sums
c
      isum=0.
      do 410 n1=1,nckt
      if(nsum(n1).eq.0)go to 410
      isum(nsum(n1))=isum(nsum(n1))+ic(n1)
410   continue
 888  format('  Time        ',16a12)
999   format(1p17e12.4)
      write(20,999)t,(vc(noc(n1))-ic(noc(n1))*
     &(ecr(noc(n1))+rps(noc(n1))),n1=1,nockt)
      write(30,999)t,(ic(noc(n1)),n1=1,nockt)
      write(40,999)t,(bz(n1),n1=1,npoint)
      write(90,999)t,(ic(noc(n1))/csac(noc(n1)),n1=1,nockt)
      write(95,999)t,(isum(n1),n1=1,nsckt)
      end if
900   continue
      end
c   ************************************************************
      subroutine lr(nckt,nsim,ntype,ncs,nsc,lc,rc,l,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   lc 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*8 (a-m,o-z)
      parameter (nc=75)
      dimension ntype(nc),ncs(nc),nsc(nc)
      dimension l(nc,nc),r(nc),lc(nc,nc),rc(nc,nc),x(nc,nc),y(nc)
      dimension v(nc),i(nc),idot(nc),vc(nc),ic(nc),icdot(nc)
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)+l(n1,n2)*idot(n2)
200   continue
      end if
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,nckt
      if(ntype(n1).eq.2)then
c
c   calculate induced voltages in open circuits
c
		vc(n1)=0.
		do 600 n2=1,nckt
		if(ntype(n2).ne.2)then
			vc(n1)=vc(n1)+lc(n1,n2)*idot(nsc(n2))
		end if
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
      end if
		icdot(n1)=idot(nsc(n1))
      end if
700   continue
      return
      end
c   ************************************************************
      subroutine build(nckt,ntype,nsim,lc,rc,l,r,x,ncs,nsc,ecr,rps)
c   ************************************************************
c
c   read icc data/cull out open circuits/invert inductance matrix
c
c   nckt is total number of 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   nsim is number of ckts after removing open ckts
c
c   lc and rc are total ind and res matricies inc'l all ckts
c
c   l and r are ind and res matricies after cullcg open ckts
c
c   rc is sum of internal ckt res (icr) and external ckt res (ecr)
c
c   x is inverse of l matrix
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   rps contains power supply equivalent resistance
c
      implicit real*8 (a-m,o-z)
      parameter (nc=75)
      dimension lc(nc,nc),rc(nc),l(nc,nc),r(nc),x(nc,nc)
      dimension ntype(nc),ncs(nc),nsc(nc)
      dimension icr(nc),ecr(nc)
      dimension rps(nc)
      rc=0.
      lc=0.
      r=0.
      l=0.
      icr=0.
      ecr=0.
      x=0.
      nsim=0
      open(1,file='cntmtx.d',status='old')
1     format(1p10e12.4)
      do 100 n=1,nckt
      read(1,1)(lc(n,n1),n1=1,nckt)
100   continue
      write(6,*)'cntmtx inductance read completed'
      read(1,1)(icr(n),n=1,nckt)
      write(6,*)'cntmtx resistance read completed'
c
c   read in external circuit resistances
c
      open(2,file='ecr.d',status='old')
      read(2,*)necr
      do 200 n=1,necr
      read(2,*)ncircuit,rexternal
      ecr(ncircuit)=rexternal
200   continue
      write(6,*)'ecr.d read completed'
c
c   write circuit inductances/resistances/time constants
c
      open(3,file='sim.p',status='old')
      do 300 n=1,nckt
      rc(n)=icr(n)+ecr(n)+rps(n)
      tau=lc(n,n)/rc(n)
3     format(i4,3e12.4)
      write(3,3)n,lc(n,n),rc(n),tau
      ncs(n)=n
      nsc(n)=n
300   continue
      l=lc
      r=rc
c
c   cull out open ciruits
c
      nopen=0
      nsim=0
      do 800 n=1,nckt
      if(ntype(n).eq.2)then
      nopen=nopen+1
c
c   eliminate row in l and r matricies
c
      do 500 nrow=nsim+1,nckt-nopen
      r(nrow)=r(nrow+1)
      do 400 ncol=1,nckt-nopen+1
      l(nrow,ncol)=l(nrow+1,ncol)
400   continue
500   continue
c
c   eliminate column in l matrix
c
      do 700 nrow=1,nckt-nopen
      do 600 ncol=nsim+1,nckt-nopen
      l(nrow,ncol)=l(nrow,ncol+1)
600   continue
700   continue
      go to 800
      end if
c
c   circuit is to be included
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
      nsim=nsim+1
      nsc(n)=nsim
      ncs(nsim)=n
800   continue
      write(6,*)'Number total, open, sim ckts =',nckt,nopen,nsim
      call invert(l,x,nsim)
      write(6,*)'matrix inversion completed'
      close(1)
      close(2)
      close(3)
      open(1,file='l.d',status='old')
      open(2,file='r.d',status='old')
      open(3,file='c.d',status='old')
      write(2,97)(r(n1),n1=1,nsim)
      do 98 n2=1,nsim
      write(1,97)(l(n1,n2),n1=1,nsim)
97    format(1p10e12.4)
98    continue
99    format(i4)
      write(3,99)(ncs(n1),n1=1,nsim)
      close(1)
      close(2)
      close(3)
      return
      end
c   ************************************************************
      subroutine fnull(nckt,ncc,ncomp,wcomp,icomp,dbdi,ic,bz,
     &npoint,bzref)
c   ************************************************************
      implicit real*8 (a-m,o-z)
      parameter(nc=75)
      parameter(nump=5)
c
c   determine required compensation currents for field null
c
c   nckt is total number of circuits
c   ncc is number of compensating circuits
c   ncomp array contains compensating circuit numbers
c   wcomp array contains weighting for each circuit
c   icomp array contains required currents
c   bzref contains reference value for Bz at null target
c   dbdi array contains Bz/amp data, first column is for null target
c   bz array contains Bz (gauss) at target + other points of interest
c
c
      dimension ncomp(10),wcomp(10),icomp(10),dbdi(nc,nump),bz(nump)
      dimension ic(nc)
c
c   calculate stray field at null target (point no. 1)
c
      sum=0.
      do 100 n=1,nckt
      sum=sum+ic(n)*dbdi(n,1)
100   continue
      bz(1)=sum*1.e4
c
c   subtract field due to compensating circuits
c
      do 200 n=1,ncc
      sum=sum-ic(ncomp(n))*dbdi(ncomp(n),1)
200   continue
c
c   determine compensating currents based on weighting
c
      do 300 n=1,ncc
      icomp(n)=(bzref-sum)*wcomp(n)/dbdi(ncomp(n),1)
300   continue
c
c   calculate field at other specified points of interest
c
      do 400 n=2,npoint
      sum=0.
      do 350 n1=1,nckt
      sum=sum+ic(n1)*dbdi(n1,n)
350   continue
      bz(n)=sum*1.e4
400   continue
      return
      end
c   ************************************************************
      subroutine invert(a,x,n)
c   ************************************************************
      implicit real*8 (a-h,o-z)
      parameter(nc=75)
c
c   compute the inverse of a matrix
c
c   reference Network Analysis Theory & Computer Methods
c   Jenkins & Watkins
c   Prentice-Hall, 1974
c   p. 86
c
c   a - input n x n matrix
c   x - inverse of a
c   n - order of a
c   eps - zero tolerance
c
      dimension a(nc,nc),x(nc,nc),c(nc,2*nc)
c
c     initialize subroutine values
c
      data eps/1.0e-12/
      n1=2*n
c
c   load a matrix and construct identity matrix
c
      do 10 i=1,n
      do 10 j=1,n
10    c(i,j)=a(i,j)
      do 30 i=1,n
      do 30 j=1,n
      if(i.eq.j)go to 20
      c(i,j+n)=0.
      go to 30
20    c(i,j+n)=1.
30    continue
c
c   invert matrix
c
      do 130 ip=1,n
c
c   find pivot element in column ip
c
      im=ip
      if(ip.ge.n)go to 50
      ist=ip+1
      do 40 i=ist,n
      if(abs(c(im,ip)).ge.abs(c(i,ip)))go to 40
      im=i
40    continue
50    if(abs(c(im,ip)).ge.eps)go to 70
c
c   near diagonal element flag
c
      write(6,60)ip,c(im,ip)
60    format(1h0,17hdiagonal element ,i2,2h =,1pe13.5)
      if(c(im,ip).eq.0.)return
70    if(im.eq.ip)go to 90
c
c   interchance rows to position pivot element
c
      do 80 j=ip,n1
      cl=c(ip,j)
      c(ip,j)=c(im,j)
80    c(im,j)=cl
c
c   find multiplication constant,set c(i,i)=1
c
90    cl=c(ip,ip)
      c(ip,ip)=1.
c
c   divide element in row by pivot element
c
      do 100 j=ist,n1
100   c(ip,j)=c(ip,j)/cl
c
c   zero column of pivot element
c
      do 120 i=1,n
      if(i.eq.ip)go to 120
      ip1=ip+1
      cl=c(i,ip)
      c(i,ip)=0.
      do 110 j=ip1,n1
110   c(i,j)=c(i,j)-cl*c(ip,j)
120   continue
130   continue
c
c   load inverse matrix into x
c
      m=n+1
      do 140 i=1,n
      do 140 j=m,n1
      k=j-n
140   x(i,k)=c(i,j)
      return
      end
