; 28, 153

FIRS_MOD = 0.5*[[1, -(1./sqrt(3.)), -(1./sqrt(3)),  (1/sqrt(3))], $
                [1, -(1./sqrt(3.)),  (1./sqrt(3)), -(1/sqrt(3))], $
                [1,  (1./sqrt(3.)), -(1./sqrt(3)), -(1/sqrt(3))], $
                [1,  (1./sqrt(3.)),  (1./sqrt(3)),  (1/sqrt(3))]] 

FIRS_MOD2 = firs_mod
firs_mod2[1:3,*]=-firs_mod2[1:3,*]

;
;  let's put in a flat field/gain error in beam 2 mod and then demod
; 

err=firs_mod*0.
percent_error=1.
for i=0,3 do begin 
   for j=0,3 do begin
      err(i,j)=firs_mod2(i,j)* randomu(seed)*percent_error/100.
   endfor
endfor

firs_mod2+=err
firs_demod=invert(firs_mod)
firs_demod2=invert(firs_mod2)

;
; Fake IQUV parameters
; 

n=100
w = findgen(n)
x=(w-n/2)/8.

i = 1.-0.7*exp( - x*x)
v=deriv(i)/3
q=-deriv(v)*.3
u=-2*q

;q*=0
; u*=0
;v*=0
;
;  realizations of time varying signals
;

nmeas=10 ; number of pol state cycles for measurement
dt=1./nmeas/4  ; one measurement state timescale
tau=dt/4 & amp=.5 ; time scale for random solar variablility in units of 1 cycle, and amplitude
f=exp(-tau/dt)*amp ; how much to change IQUV between states
;
; dual or single beam setup
;
wt=[1.0,.0]                     ; single beam
wt=[1.0,1.]                    ; dual beam
wt/=total(wt)

outvec=fltarr(n,nmeas,4)
for j=0,n-1 do begin ; wavelength index
   invec=[i[j],q[j],u[j],v[j]]

   for k=0,nmeas-1 do begin  ; number of cycles
      state=fltarr(4)
      state2=fltarr(4)
      dstate=fltarr(4)
      dstate2=fltarr(4)
      for m=0,3 do begin  ; number of states per comlete pol cycle
;
; for seeing induced changes, use only one random num for all IQUV
;
         di=randomn(seed)*f*i[j]
         dq=randomn(seed1)*f*q[j]
         du=randomn(seed2)*f*u[j]
         dv=randomn(seed3)*f*v[j]
;
;  secular changes
;
         di += 0.5*k/nmeas*i[j]  ; increasing with time
         dq += 0.5*k/nmeas*q[j]  ; increasing with time
         du += 0.5*k/nmeas*u[j]  ; increasing with time
         dv += 0.5*k/nmeas*v[j]  ; increasing with time
; vector
         vec=[di,dq,du,dv]      ; build noisy vector at each state
; beam 1
         dstate[m]= total(firs_mod(*,0)*vec)
; beam 2
         dstate2[m]= total(firs_mod2(*,0)*vec)
      endfor
      state=firs_mod##invec + dstate
      state2=firs_mod2##invec + dstate2
      
      outvec[j,k,*]=wt[0]*firs_demod## state + wt[1]*firs_demod2##state2
      if(j eq -1) then begin
         print
         print,'state',transpose(state)-dstate
         print,'dstate',dstate
         print,'demodulated '
         print,firs_demod ##state
         print
         stop
      end
   endfor
endfor

!p.multi=[0,2,2]

for kk=0,3 do begin 
   y=mean(outvec[*,*,kk],dim=1)
   plot,w,y,lines=3,ps=-1
   case kk of
      0:z=i
      1:z=q
      2:z=u
      3:z=v
   endcase
   oplot,w,z,lines=0
   s=stdev(y-z)
   xyouts,10,!y.crange(0)+0.1*(!y.crange(1)-!y.crange(0)),string(s,form='(f8.4)')
endfor

!p.multi=0



end