dblGauss

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  function vis = doubleGauss(d, nu, par, beam, pixSize)


  this function should produce the visibilities of a double
  gaussian.

  d has your u,v,w information, par is a vector with the following
  information:
  par = [Amp_30GHz, alpha, theta, l_centroid, m_centroid, sig_l, sig_m]

  where theta is the angle of inclination of the axes of the gaussian.

  beam is a vector, [1 by 16], with the beam value at the centroid
  l,m for all frequencies.

  pixSize is necessary for proper normalization (in radians)


  we solve for this analytically, to see my long and drawn our
  solution, look in my lab book on pages 72-77.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function vis = dblGauss(d, nu, par, beam, pixSize )
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function vis = doubleGauss(d, nu, par, beam, pixSize)
0006 %
0007 %
0008 %  this function should produce the visibilities of a double
0009 %  gaussian.
0010 %
0011 %  d has your u,v,w information, par is a vector with the following
0012 %  information:
0013 %  par = [Amp_30GHz, alpha, theta, l_centroid, m_centroid, sig_l, sig_m]
0014 %
0015 %  where theta is the angle of inclination of the axes of the gaussian.
0016 %
0017 %  beam is a vector, [1 by 16], with the beam value at the centroid
0018 %  l,m for all frequencies.
0019 %
0020 %  pixSize is necessary for proper normalization (in radians)
0021 %
0022 %
0023 %  we solve for this analytically, to see my long and drawn our
0024 %  solution, look in my lab book on pages 72-77.
0025 %
0026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0027 
0028 % first we get our values and parse them.
0029 r = sqrt(par(4).^2 + par(5).^2);  % how far from beam center the
0030                                   % beast is  (in degrees)
0031 amp   = par(1);
0032 alpha = par(2);
0033 theta = par(3)*pi/180;
0034 l_cen = par(4)*pi/180;
0035 m_cen = par(5)*pi/180;
0036 sig_l = par(6)/60*pi/180;
0037 sig_m = par(7)/60*pi/180;
0038 
0039 cosTheta = cos(theta);
0040 sinTheta = sin(theta);
0041 
0042 l_cen_prime = l_cen*cosTheta - m_cen*sinTheta;
0043 m_cen_prime = l_cen*sinTheta + m_cen*cosTheta;
0044 
0045 % since our beam is symmetric, it doesn't matter whether it's in
0046 % regular or primed coordinates
0047 A = amp*(nu/30e9).^(-alpha);
0048 %norm = 2*pi*sig_l*sig_m/(pixSize)^2;
0049 %A = A*norm;
0050 if(isstruct(beam))
0051   f = find(abs(beam.r-r) == min(abs(beam.r - r)));
0052   bmVal = beam.beam(f,:);
0053 else
0054   bmVal = beam;
0055 end
0056 
0057 %bmVal = 1;
0058 
0059 % now we calculate the constant factors
0060 F1 = A.*bmVal;
0061 
0062 F2 = sig_l^2.*(d.u{1}*cosTheta - d.v{1}*sinTheta).^2 + sig_m.^2.*(d.u{1}*sinTheta + d.v{1}*cosTheta).^2;
0063 F2 = exp(-2*pi^2*F2);
0064 
0065 F3 = -2*pi*[(d.u{1}*cosTheta - d.v{1}*sinTheta)*l_cen_prime + (d.u{1}*sinTheta + d.v{1}*cosTheta)*m_cen_prime];
0066 
0067 if(ndims(F2)==3)
0068   F1 = reshape(F1, [1 1 16]);
0069   F1 = repmat(F1, [size(F2,1) size(F2,2) 1]);
0070 elseif(ndims(F2)==2)
0071   F1 = repmat(F1, [size(F2,1) 1]);
0072 end
0073 F = F1.*F2;
0074 vis{1} = complex(F.*cos(F3), F.*sin(F3));
0075 
0076 return;
0077 
0078 % par = [1 1 0 0 0 1 2]
0079 % beam = ones(1,16)
0080 
0081 %par = mosaic_srcfit('dblGauss', peak, [0 2000], ad, ads, o, p,beam);