; 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