mosaic_srcfit

 [par,p, loc]=mosaic_srcfit(type,peak,uvcut,ad,ads,o,p,beam,g)

 Do a model fit to multi pointing visibility data

 type = 'cl' or 'ps' for clusters or point source
 peak = location and peak flux of source
 uvcut = uv range to fit to
 ad/ads = image definition structures
 o = observation data (only o.nu used)
 p = pointing data (for px/y,u,v,vis,var)
 beam = primary beam at each band

 par = model parameters
 p = multi pointing data with model subtracted (from all pointing and
 uvr)
 loc = location in J2000 string coordinates.

0001 function [par,p, loc]=mosaic_srcfit(type,peak,uvcut,ad,ads,o,p,beam)
0002 % [par,p, loc]=mosaic_srcfit(type,peak,uvcut,ad,ads,o,p,beam,g)
0003 %
0004 % Do a model fit to multi pointing visibility data
0005 %
0006 % type = 'cl' or 'ps' for clusters or point source
0007 % peak = location and peak flux of source
0008 % uvcut = uv range to fit to
0009 % ad/ads = image definition structures
0010 % o = observation data (only o.nu used)
0011 % p = pointing data (for px/y,u,v,vis,var)
0012 % beam = primary beam at each band
0013 %
0014 % par = model parameters
0015 % p = multi pointing data with model subtracted (from all pointing and
0016 % uvr)
0017 % loc = location in J2000 string coordinates.
0018 
0019 
0020 % Reduce input 2D beam to 1D (assume axi-symmetric)
0021 % if we use perfect beam we won't have
0022 % to worry about orientation.
0023 n=floor(size(beam,1)/2)+1;
0024 bm.beam=squeeze(beam(n:end,n,:));
0025 bm.r=ads.t_val_deg(n:end)';
0026 
0027 % find x and y location for each pointing
0028 xp=peak.x-p.px;
0029 yp=peak.y-p.py;
0030 
0031 % Find close pointings
0032 r=pyth(xp,yp);
0033 close=find(r<0.15);
0034 if(isempty(close))
0035     close = find(r == min(r));
0036     display('no ultra close pointings');
0037     display(sprintf('closest pointing at: %f degrees', min(r)));
0038     if(min(r)>.15)
0039         % closest pointing is where primary beam is less than 25%.
0040         display('really really really far');
0041         display('can not fit');
0042         par = zeros(1,4);
0043         p = p;
0044         loc = [];
0045         return;
0046     end
0047 end
0048 near = find(r<(ads.Field_size_deg/2));
0049 disp(sprintf('%d near pointings', length(near)));
0050 
0051 if(~isfield(p, 'newFit'))
0052   p.newFit = zeros(1,length(p.n));
0053 end
0054 p.newFit(near) = 1;
0055 
0056 
0057 % find matching uvcut range
0058 if(close)
0059 %%% Begin New Section
0060 % note change in all stored uvrange maps
0061   rm=[];
0062   for i=1:length(p.mos)
0063     % if remap not empty, then map is not up-to-date
0064     % and should be removed
0065     if (isempty(p.mos{i}.remap))
0066       p.mos{i}.remap=[1:10];
0067 %      p.mos{i}.remap=close;
0068     else
0069       rm=[rm i];
0070     end
0071   end
0072   p.mos(rm)=[];
0073 %%% End New Section
0074  disp(sprintf('%d close pointings',length(close)));
0075 else
0076   disp('Cannot fit out source - no close pointings')
0077   par=[];
0078   return
0079 end
0080 
0081 % cut down to close pointings
0082 pc=rmfield(p,'mos');
0083 if(length(close)>1)
0084     pc=structcut(pc,close);
0085 end
0086 pn = rmfield(p, 'mos');
0087 if(length(near)>1)
0088     pn = structcut(pn, near);
0089 end
0090 switch type
0091   % Point Source fit
0092  case 'ps'
0093   par=lsqnonlin(@(x) mchisq(pc,psvisfunc(x,o.nu,pc,bm),uvcut), ...
0094         [peak.flux,0,peak.x,peak.y]);
0095   mod=psvisfunc(par,o.nu,pn,bm);
0096   
0097   % Point source with fixed spectral index
0098  case 'psAlpha'
0099   par=lsqnonlin(@(x) mchisq(pc,psvisfunc(x,o.nu,pc,bm, peak.vla.vlaF),uvcut), ...
0100         [peak.flux,NaN,peak.x,peak.y]);  
0101   mod=psvisfunc(par,o.nu,pn,bm, peak.vla.vlaF);    
0102     
0103   % Cluster with SZ spectral shape
0104  case 'cl'
0105   % Fit 1 - fixed width - Starting width val is 1 arcmin in radians
0106   par=[peak.flux,3e-4,peak.x,peak.y];
0107   lo=[-Inf,3e-4-1e-9,-Inf,-Inf]; up=[+Inf,3e-4+1e-9,+Inf,+Inf];
0108   
0109   par=lsqnonlin(@(x) mchisq(pc,clvisfunc(x,o.nu,pc,bm),uvcut),par,lo,up);
0110   
0111   % Fit 2 - free width
0112   par=lsqnonlin(@(x) mchisq(pc,clvisfunc(x,o.nu,pc,bm),uvcut),par);
0113   
0114   % calc best fit model for all pointings and basline lengths
0115   mod=clvisfunc(par,o.nu,p,bm);
0116   
0117   
0118   % Elliptical Cluster with SZ spectral shape
0119   case 'cle'
0120     % Fit 1 - fixed width - Starting width val is 1 arcmin in radians
0121     % also fixed is ratio of elliptical axes.  Assume circular
0122     par=[peak.flux,3e-4,peak.x,peak.y,1,0];
0123     lo=[-Inf,3e-4-1e-9,-Inf,-Inf,1-1e-9,-1e-9]; 
0124     up=[+Inf,3e-4+1e-9,+Inf,+Inf,1+1e-9,+1e-9];
0125     tic
0126     par=lsqnonlin(@(x) mchisq(pc,clevisfunc(x,o.nu,pc,bm),uvcut),par,lo,up);
0127 
0128     % Fit 2 - free width, axes ratio, and ellipse orientation
0129     par=lsqnonlin(@(x) mchisq(pc,clevisfunc(x,o.nu,pc,bm),uvcut),par);
0130     
0131     % calc best fit model for all pointings and basline lengths
0132     mod=clvisfunc(par,o.nu,p,bm);
0133     toc   
0134     
0135   % "Cluster" with power law spectra
0136   case 'cl2'
0137     % Fit 1 - fixed width - Starting width val is 1 arcmin in radians
0138     par=[peak.flux,0,3e-4,peak.x,peak.y];
0139     lo=[-Inf,-Inf,3e-4-1e-9,-Inf,-Inf]; up=[+Inf,+Inf,3e-4+1e-9,+Inf,+Inf];
0140     
0141     par=lsqnonlin(@(x) mchisq(pc,cl2visfunc(x,o.nu,pc,bm),uvcut),par,lo,up);
0142 
0143     % Fit 2 - free width
0144     %par=lsqnonlin(@(x) mchisq(pc,cl2visfunc(x,o.nu,pc,bm),uvcut),par);
0145     
0146     % calc best fit model for all pointings and basline lengths
0147     mod=cl2visfunc(par,o.nu,pn,bm);
0148 
0149 
0150   % "Double Gaussian"
0151  case 'dblGauss'
0152    OPTIONS = optimset('Display', 'off', 'MaxIter', 200000, ...
0153              'MaxFunEvals', 20000, 'TolX', 0.00005, 'TolFun', 0.00005);
0154   par = [peak.flux, -0, 0, peak.x, peak.y, 0.1,0.1];
0155   if(peak.flux>0)
0156       lo  = [0, -2, -180, -3, -3, 0, 0];
0157       hi  = [Inf, 2, 180, 3, 3, 180, 180];
0158   else
0159       lo  = [-Inf, -2, -180, -3, -3, 0, 0];
0160       hi  = [0, 2, 180, 3, 3, 180, 180];
0161   end
0162   
0163   if(isfield(peak, 'vla'))
0164     % this means we're using priors from NVSS
0165     [par, resnorm, residual, exitflag, output] = lsqnonlin(@(x) mchisq(pc, dblGaussVisFunc(pc,o.nu,x, bm,peak), uvcut),par,lo,hi,OPTIONS);
0166     par(2) = -log(peak.vla.vlaF/(1000*par(1)))/(log(1.4e9/mean(o.nu)));
0167   else
0168     [par, resnorm, residual, exitflag, output] = lsqnonlin(@(x) mchisq(pc, dblGaussVisFunc(pc,o.nu,x, bm), uvcut),par,lo,hi,OPTIONS);
0169   end
0170   mod=dblGaussVisFunc(pn,o.nu,par,bm);
0171 end
0172 
0173 
0174 % Subtract best fit model
0175 for i=1:length(pn.px)
0176   modn = mod{i};
0177   modn(isnan(modn)) = 0;
0178   pn.vis{i}=pn.vis{i}-modn;
0179 end
0180 
0181 for m=1:length(near)
0182   p.vis{near(m)} = pn.vis{m};
0183 end
0184 
0185 % Put in the best fit parameters in coordinates that are useful to us.
0186 
0187 if(length(par)==7)
0188   [RA DEC] = fieldoff_to_radec(par(4), par(5), o.rafc*15, o.decfc);
0189 else
0190   [RA DEC] = fieldoff_to_radec(par(3), par(4), o.rafc*15, o.decfc);
0191 end
0192 [loc.RA loc.DEC] = radec2str(RA/15, DEC);
0193 
0194 if(isnan(par(2)))
0195     mnu = mean(o.nu);
0196     alpha = log(peak.vla.vlaF/par(1))/(log(1.4e9/mnu));
0197     par(2) = alpha;
0198 end
0199 
0200 return
0201 
0202 
0203 
0204 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0205 function V=clevisfunc(par,nu,p,bm)
0206 % Fit visibility data to exponential profile
0207 % Think this is FT of beta model with beta=4/3 as favored by Holder et al
0208 %
0209 % par(1) = Visibility amp at zero spacing and band center -
0210 % positive as SZ spectral shape includes the negative
0211 % par(2) = Exp width
0212 % par(3:4) = source location on field
0213 % par(5) = a/b (ratio of major/minor axes)
0214 % par(6) = angle [deg]
0215 
0216 % Calc relative SZ flux at each freq
0217 % Normalize to unity at band center
0218 nu=[nu,mean(nu)];
0219 Tcmb=2.726;
0220 fx=sz_fx(nu); delTsz=Tcmb.*fx;
0221 dIdT=planck_dIdT(nu,Tcmb);
0222 delIsz=delTsz.*dIdT;
0223 delIsz=delIsz./delIsz(end);
0224 d=delIsz(1:end-1);
0225 % Combine with fit amp - the zero spacing flux
0226 s=cvec(par(1).*d);
0227 
0228 % get dot products and beam values
0229 [ap,bp]=vissubfunc(par(3),par(4),p,bm);
0230 
0231 % convert deg to rad
0232 rot=par(6)*pi/180;
0233 
0234 % for each pointing
0235 for i=1:length(p.px)
0236   % calculate ellpise ring
0237   x=p.u{i}*cos(rot)+p.v{i}*sin(rot);
0238   y=-p.u{i}*sin(rot)+p.v{i}*cos(rot);
0239   [theta,r]=cart2pol(x,y);
0240 
0241   % expand out flux value
0242   sp=s(p.band{i});
0243 
0244   % Compute vis amp rel to zero spacing for this baseline length
0245   vp=exp(-2.*pi.*elip(theta,par(2),par(5)).*r);
0246 
0247   % Expected visibility
0248   V{i}=sp.*vp.*bp{i}.*ap{i};
0249 end
0250 
0251 return
0252 
0253 
0254 %%%%%%%%%%%%%%%%%%%%%%%
0255 function e=elip(theta,g,s)
0256 
0257 % elliptical scaling of exponential width
0258 %
0259 % exp(2*pi*f(theta)*r)
0260 %
0261 % f(theta)=g*(ro/r1)
0262 %   where ro is defined to be major axis
0263 %   and r1 is the radius of any point on the reference ellipse defined
0264 %   by (ax,by) relative to the origin
0265 %
0266 % ellipse function defined by
0267 %
0268 % u^2/ax^2+v^2/by^2=1
0269 
0270 % to keep scaling correct ro must be reference to major axis
0271 % if ax is larger, then we are ok because theta is calculated relative
0272 % to x.  if by is larger, then we must subtract by pi/2 to make angle
0273 % relative to positive y axis
0274 if (s>=1)
0275   a=s;
0276   b=1;
0277 else
0278   a=1;
0279   b=s;
0280   theta=pi/2-theta;
0281 end
0282   
0283 arg1=1/a^2+(sin(theta).^2./((b^2)*(cos(theta).^2)));
0284 arg2=1-(b^2)/(a^2);
0285 
0286 % r1/ro
0287 e=(1/a)*sqrt((1./arg1)*arg2+b^2);
0288 
0289 % need to inverti
0290 e=(1./e)*g;
0291 
0292 
0293 
0294 
0295 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0296 function V=clvisfunc(par,nu,p,bm)
0297 % Fit visibility data to exponential profile
0298 % Think this is FT of beta model with beta=4/3 as favored by Holder et al
0299 %
0300 % par(1) = Visibility amp at zero spacing and band center -
0301 % positive as SZ spectral shape includes the negative
0302 % par(2) = Exp width
0303 % par(3:4) = source location on field
0304 
0305 % Calc relative SZ flux at each freq
0306 % Normalize to unity at band center
0307 nu=[nu,mean(nu)];
0308 Tcmb=2.726;
0309 fx=sz_fx(nu); delTsz=Tcmb.*fx;
0310 dIdT=planck_dIdT(nu,Tcmb);
0311 delIsz=delTsz.*dIdT;
0312 delIsz=delIsz./delIsz(end);
0313 d=delIsz(1:end-1);
0314 % Combine with fit amp - the zero spacing flux
0315 s=cvec(par(1).*d);
0316 
0317 % get dot products and beam values
0318 [ap,bp]=vissubfunc(par(3),par(4),p,bm);
0319 
0320 % for each pointing
0321 for i=1:length(p.px)
0322 
0323   % expand out flux value
0324   sp=s(p.band{i});
0325 
0326   % Compute vis amp rel to zero spacing for this baseline length
0327   vp=exp(-2*pi*par(2)*p.uvr{i});
0328 
0329   % Expected visibility
0330   V{i}=sp.*vp.*bp{i}.*ap{i};
0331 end
0332 
0333 return
0334 
0335 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0336 function V=cl2visfunc(par,nu,p,bm)
0337 % Fit visibility data to exponential profile with pow law spec index
0338 % par(1) = Visibility amp at zero spacing and band center -
0339 % positive as SZ spectral shape includes the negative
0340 % par(2) = spectral index
0341 % par(3) = Exp width
0342 % par(4:5) = source location on field
0343 
0344 % Flux is flux at mean freq
0345 mnu=mean(nu);
0346 
0347 % Calc src flux for each freq
0348 s=cvec(par(1).*(nu/mnu).^par(2));
0349 
0350 % get dot products and beam values
0351 [ap,bp]=vissubfunc(par(4),par(5),p,bm);
0352 
0353 % for each pointing
0354 for i=1:length(p.px)
0355 
0356   % expand out flux value
0357   sp=s(p.band{i});
0358 
0359   % Compute vis amp rel to zero spacing for this baseline length
0360   vp=exp(-2*pi*par(3)*p.uvr{i});
0361 
0362   % Expected visibility
0363   V{i}=sp.*vp.*bp{i}.*ap{i};
0364 end
0365 
0366 return
0367 
0368 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0369 function V=psvisfunc(par,nu,p,bm,Svla)
0370 % Expected vis for a point source
0371 %
0372 % par(1) = flux
0373 % par(2) = spectral index
0374 % par(3:4) = source location on field in deg
0375 % Svla is optional, it's the flux from the VLA.
0376 % Flux is flux at mean freq
0377 mnu=mean(nu);
0378 
0379 if(exist('Svla'))
0380   % we're fixing a spectral index from SZA/VLA
0381   alpha = -log(Svla/(1000*par(1)))/(log(1.4e9/mnu));
0382 else
0383   alpha = par(2);
0384 end
0385 
0386 % Calc src flux for each freq
0387 s=cvec(par(1).*(nu/mnu).^(-alpha));
0388 
0389 % get dot products and beam values
0390 [ap,bp]=vissubfunc(par(3),par(4),p,bm, nu);
0391 
0392 % for each pointing
0393 for i=1:length(p.px)
0394   
0395   % expand out flux value
0396   sp=s(p.band{i});
0397   
0398   % Expected visibility is src flux * beam taper * sin/cos dot product
0399   % of src position and baseline dir
0400   V{i}=sp.*bp{i}.*ap{i};  
0401 end
0402 
0403 return
0404 
0405 
0406 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0407 function V =dblGaussVisFunc(p,nu,par,bm, peak)
0408 %  Fit visibility data to double Gaussian with power law spec index
0409 %  par(1) = Visibility amp at zero spacing and band center -
0410 %  par(2) = spectral index
0411 %  par(3) = angle of orientation of ellipse
0412 %  par(4,5) = location of center
0413 %  par(6,7) = sigma in each direction.
0414 
0415 % calc src location in each pointing
0416 px=-(par(4)-p.px); py=-(par(5)-p.py);
0417 pr=pyth(px,py);
0418 
0419 % for each band calc beam values at this location for each pointing
0420 for i=1:size(bm.beam,2)
0421   b(:,i)=interp1(bm.r,bm.beam(:,i),pr);
0422 end
0423 % no use for theoretical beam...interp is fast.
0424 
0425 mnu = mean(nu);
0426 
0427 for m=1:length(p.px)
0428   % calculate parameter for this field
0429   thisPar = par;
0430   if(exist('peak'))
0431     thisPar(2) =  -log(peak.vla.vlaF/(1000*par(1)))/(log(1.4e9/mnu));
0432     if(isfield(peak.vla, 'vlaAngle'))
0433       thisPar(3) =  peak.vla.vlaAngle;
0434     end
0435     if(isfield(peak.vla, 'vlaRatio'))
0436       thisPar(6) = peak.vla.vlaRatio*thisPar(7);
0437     end
0438   end
0439   thisPar(4) = px(m);
0440   thisPar(5) = py(m);
0441   d.u{1} = p.u{m};
0442   d.v{1} = p.v{m};
0443   vis  = dblGauss(d, nu, thisPar, b(m,:));
0444   V{m} = vis{1};
0445 end
0446 
0447 return;
0448 
0449 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0450 function [ap,bp]=vissubfunc(x,y,p,bm, nu)
0451 
0452 % calc src location in each pointing
0453 px=x-p.px; py=y-p.py;
0454 pr=pyth(px,py);
0455 
0456 % for each band calc beam values at this location for each pointing
0457 for i=1:size(bm.beam,2)
0458   b(:,i)=interp1(bm.r,bm.beam(:,i),pr);
0459 end
0460 % no use for theoretical beam...interp is fast.
0461 
0462 % convert to rad
0463 px=px*pi/180;
0464 py=py*pi/180;
0465 
0466 % for each pointing
0467 for i=1:length(px)
0468   % Dot product of src location and baseline vectors for this pointing
0469   n=sqrt(1-px(i)^2-py(i)^2);
0470   x=2*pi*(px(i)*p.u{i}+py(i)*p.v{i}+p.w{i}*(n-1));
0471   ap{i}=complex(cos(x),sin(x));
0472   
0473   % expand out beam values
0474   bb=b(i,:)/n;
0475   bp{i}=bb(p.band{i});
0476 end
0477 
0478 return
0479 
0480 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0481 function x=mchisq(p,modvis,uvcut)
0482 % x=mchisq(p,modvis,uvcut)
0483 %
0484 % Make (data-model)/error for lsqnonlin
0485 
0486 disp('mchisq')
0487 
0488 % concat pointings
0489 vis=vertcat(p.vis{:});
0490 err=sqrt(vertcat(p.var{:}));
0491 uvr=vertcat(p.uvr{:});
0492 modvis=vertcat(modvis{:});
0493 
0494 % vectorize
0495 vis=cvec(vis);
0496 uvr=cvec(uvr);
0497 modvis=cvec(modvis);
0498 err=cvec(err);
0499 
0500 % apply uvr cut
0501 ind=uvr>uvcut(1)&uvr<uvcut(2);
0502 vis=vis(ind);
0503 modvis=modvis(ind);
0504 err=err(ind);
0505 
0506 % remove flagged data
0507 i=~isnan(vis);
0508 vis=vis(i);
0509 modvis=modvis(i);
0510 err=err(i);
0511 
0512 % delta obs/model
0513 D=vis-modvis;
0514 
0515 D=[real(D);imag(D)];
0516 err=[err;err];
0517 
0518 
0519 x=D./err;
0520 
0521 % remove NaN or Inf data
0522 x(find(isnan(x))) = [];
0523 x(find(isinf(x))) = [];
0524 
0525 return