function ssccd, sysc, Ts , tau

;  sysd = ssccd(sysc,tau) computes a state realization for systems
;  with CCD-based measurements based on the continuous time system
;  sysc (sysc.A,sysc.B,sysc.C,sysc.D) for sampling time Ts and delay tau
;
;  The output is a discrete-time state variable model with structure
;	sysd.a, sysd.b, sysd.c, sysd.d.

  compile_opt idl2
  IF tau EQ 0 THEN tau = Ts*1.e-4 ; small delay so that the program runs
  
  nda = N_TAGS(sysc)
  case nda of
     1: na = 0
     4: begin
        sza = SIZE(sysc.a)
        if sza[0] eq 0 then na = 1 else na = sza[1]
        a = sysc.a
        b = sysc.b
        c = sysc.c
     end
  endcase

  szd = SIZE(sysc.d)
  d = sysc.d

  case szd[0] of
     0: begin
        mb = 1
        pc = 1
     end
     1: begin
        mb = szd[1]
        pc = 1
     end
     2: begin
        mb = szd[1]
        pc = szd[2]
     end
  endcase

  if na eq 0 then begin
     abig = [ [ FLTARR(pc+mb,pc) ] , [ FLTARR(pc+mb,mb)] ]
  endif else begin
     abig = [ [ FLTARR(pc,pc) , c , FLTARR(mb,pc)] , [ FLTARR(pc,na) , a , b ] , [ FLTARR(pc+na+mb,mb) ] ]
  endelse

  ind1 = INDGEN(pc)
  ind1b = pc-1
  if na gt 0 then begin
     ind2a = pc
     ind2b = pc+na-1
     ind2 = INDGEN(na)+pc
  endif
  ind3 = INDGEN(mb)+pc+na
  ind3a = pc+na
  ind3b = pc+na+mb-1

  fTs = EXPM(abig*Ts)           ;

  if na gt 0 then begin
     ad = fTs[ind2a:ind2b,ind2a:ind2b]
     bd = fTs[ind3a:ind3b,ind2a:ind2b]
  endif

  if tau ne 0 then begin
     el = CEIL(tau/Ts)
     beta = el*Ts-tau
     fbeta = EXPM(abig*beta)
     if na gt 0 then begin
        cccd0 = fbeta[ind2a:ind2b,0:ind1b] ;
        cccd1 = fTs[ind2a:ind2b,0:ind1b]-fbeta[ind2a:ind2b,0:ind1b]
     endif
     dccd0 = fbeta[ind3a:ind3b,0:ind1b]+beta*d
     dccd1 = fTs[ind3a:ind3b,0:ind1b]-fbeta[ind3a:ind3b,0:ind1b]+(Ts-beta)*d
     if na gt 0 then begin
        br1a = [ [[ad],[cccd1]],FLTARR((el+1)*pc,pc+na) ]
        case el of
           2: begin
              br2a = [ [[cccd0],[FLTARR(na,(el-1)*pc)]],IDENTITY(el*pc),FLTARR(pc,el*pc)]
              bccd = [[bd],[dccd1],[dccd0],[FLTARR(mb,(el-1)*pc)]]
           end
           1: begin
              br2a = [ [cccd0],IDENTITY(el*pc),FLTARR(pc,el*pc)]
              bccd = [[bd],[dccd1],[dccd0]]
           end
           0: begin
              br2a = [ cccd0 ]
              bccd = [[bd],[dccd1],[dccd0]]
           end
        endcase
        accd = [[br1a],[br2a]]
        cccd = [FLTARR(el*pc+na,pc), IDENTITY(pc)]
        dccd = FLTARR(mb, pc)
     endif else begin
        br1a = FLTARR((el+1)*pc,pc+na)
        case el of
           2: begin
              br2a = [IDENTITY(el*pc),FLTARR(pc,el*pc)]
              accd = [[br1a],[br2a]]
              bccd = [[dccd1],[dccd0],[FLTARR((el-1)*pc,mb)]]
           end
           1: begin
              br2a = [IDENTITY(el*pc),FLTARR(pc,el*pc)]
              accd = [[br1a],[br2a]]
              bccd = [[dccd1],[dccd0]]
           end
           0: begin
              accd = br1a
              bccd = [[dccd1],[dccd0]]
           end
        endcase
        cccd = [FLTARR(el*pc+na,pc), IDENTITY(pc)]
        dccd = FLTARR(pc,mb)
     endelse
  endif else begin
     dccd1 = fTs[ind3,ind1]+Ts*d ;
     if na gt 0 then begin
        cccd1 = fTs[ind2a:ind2b,ind1]
        accd = [[ad,FLTARR(pc,na)],[cccd1, FLTARR(pc,pc)]]
        bccd = [[bd],[dccd1]]
        cccd = [FLTARR(na,pc),IDENTITY(pc)]
     endif else begin
        accd = FLTARR(pc,pc)
        bccd = dccd1
        cccd = [IDENTITY(pc)]
     endelse
     dccd = FLTARR(pc,mb)       ;
  endelse

  sysd = CREATE_STRUCT('a',accd,'b',bccd,'c',cccd,'d',dccd)
  return, sysd

end