polcal2

 function to generate calibration matrix for IQU data in form required by
 descart

0001 % function to generate calibration matrix for IQU data in form required by
0002 % descart
0003 function [calmat fval hessian] = polcal2(data, model)
0004 
0005 % data are corresponding values of I1,I2,Q1,Q2,U1,U2 and phi (parallactic angle)
0006 % as an n-by-7 matrix
0007 % model is a vector, [I (intensity) p (polarization fraction) theta (position angle)]
0008 % such that the model Stokes parameters are given by
0009 %   I = I
0010 %   Q = I * p * cos(2(theta + phi))
0011 %   U = I * p * sin(2(theta + phi))
0012 % calmat is the 3-by-6 Mueller matrix that is the
0013 % telescope response, i.e. if the true sky Stokes vector is S and the
0014 % observed vector is S', S' = SM
0015 % so what we want to minimise is
0016 %      chisq = sum(mod-squared(SM - S'))
0017 % as a function of M
0018 
0019 %(I,Q,U)(110000   = (I1,I2,Q1,Q2,U1,U2)
0020 %        001100
0021 %        000011)
0022 
0023 calmat0 = [1,1,0,0,0,0;...
0024            0,0,1,1,0,0;...
0025            0,0,0,0,1,1];      %starting value
0026 
0027 % call fminunc using the nested function method
0028 % fminunc treats calmat as a vector (col1,col2,col3,...col6)
0029 % Resulting hessianhas outputs ordered (I_I1 Q_I1 U_I1 I_I2 Q_I2 U_I2....)
0030 options = optimset('PlotFcns',@optimplotfval,'TolFun',1e-18);
0031 [calmat, fval, hessian] = fminunc(@chisq, calmat0,options)
0032 
0033 %[par,err,gof,stat] = matmin('chisq',calmat0,[],func,y,...)
0034 % Try lsqnonlin
0035 %[calmat,fval,hessian]=lsqnonlin(@chi,calmat0)
0036 
0037 
0038 %Now plot some checks
0039 figure
0040 k=0
0041 
0042 for z=1:3
0043     for w=1:6
0044         calmat_temp=calmat;
0045         k=k+1;
0046         %p=calmat(z,w);
0047         m=0;
0048         for l=-0.5:0.01:0.5
0049             m=m+1;
0050             calmat_temp(z,w)=calmat(z,w)+l;
0051             chisq_temp(m)=chisq(calmat_temp);
0052             x(m)=calmat(z,w)+l;
0053         end
0054         subplot(3,6,k)
0055         %keyboard
0056         plot(x,chisq_temp)
0057         hold all
0058         plot(calmat(z,w),chisq(calmat),'or')
0059         % Calculate the curvature
0060         [minval minpos]=min(chisq_temp);
0061 %        curvature(z,w)=(chisq_temp(minpos+1) + chisq_temp(minpos-1) -2*chisq_temp(minpos) )/(x(minpos+1)-x(minpos))^2;
0062         
0063     end
0064    
0065 end
0066 
0067 %curvature
0068 %  err
0069 %  gof
0070 %  stat
0071             
0072             
0073             
0074 % Use simulanneal to minimise
0075 %
0076 % options = saoptimset('PlotFcns',{@saplotbestx,...
0077 %                @saplotbestf,@saplotx,@saplotf});
0078 % [calmat,fval]=simulannealbnd(@chisq,calmat0,[],[],options)
0079 % Try lsqnonlin
0080 % [calmat2,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@chi,calmat0)
0081 
0082     function f = chisq(calmat)
0083         f = 0;
0084         npoints = length(data(:,1));
0085         for i = 1:npoints;
0086             phi = -data(i,7);
0087             modelvec = [model(1), model(1)*model(2)*cosd(2*(model(3)+phi)), model(1)*model(2)*sind(2*(model(3)+phi))]; % model IQU in telescope coords
0088             datavec = [data(i,1), data(i,2), data(i,3),data(i,4), data(i,5), data(i,6)]; %I1,I2,Q1,Q2,U1,U2
0089             chi = modelvec*calmat - datavec;
0090             %chi = calmat'*modelvec' - datavec';
0091             f = f + norm(chi)^2;
0092         end
0093     end
0094 
0095 function g = chi(calmat)
0096         g = [];
0097         npoints = length(data(:,1));
0098         for i = 1:npoints;
0099             phi = -data(i,7);
0100             modelvec = [model(1), model(1)*model(2)*cosd(2*(model(3)+phi)), model(1)*model(2)*sind(2*(model(3)+phi))]; % model IQU in telescope coords
0101             datavec = [data(i,1), data(i,2), data(i,3),data(i,4), data(i,5), data(i,6)]; %I1,I2,Q1,Q2,U1,U2
0102             
0103             chi = modelvec*calmat - datavec;
0104             g = [g,chi];
0105         end
0106 end
0107 
0108 function h = check(data,model)
0109         h = 0;
0110         npoints = length(data(:,1));
0111         for i = 1:npoints;
0112             phi = -data(i,7);
0113             modelvec = [model(1), model(1)*model(2)*cosd(2*(model(3)+phi)), model(1)*model(2)*sind(2*(model(3)+phi))]; % model IQU in telescope coords
0114             datavec = [data(i,1), data(i,2), data(i,3),data(i,4), data(i,5), data(i,6)]; %I1,I2,Q1,Q2,U1,U2
0115             h = modelvec*calmat - datavec;
0116             %chi = calmat'*modelvec' - datavec';
0117             % f= f + norm(chi)^2;
0118         end
0119     end
0120 
0121 
0122 
0123 end
0124 %return calmat