work2D

frank/Examples/AS209_1mm_FrankFit.py

@@ -0,0 +1,44 @@
+import frank
+import numpy as np 
+import matplotlib.pyplot as plt
+
+from frank.geometry import SourceGeometry
+from frank.radial_fitters import FrankFitter, FourierBesselFitter
+
+# Huang 2018 
+inc = 34.97
+pa = 85.76
+dra = 1.9e-3
+ddec = -2.5e-3
+r_out = 1.9
+
+# Frank Parameters
+n_pts = 300
+alpha = 1.3
+w_smooth = 1e-1
+
+# UVtable AS209 at 1mm with removed visibilities between .
+dir = "./"
+data_file = dir + 'AS209_continuum_prom_1chan_30s_keepflagsFalse_removed1.txt'
+
+# Loading data
+u, v, Re, Im, Weights = np.loadtxt(data_file, unpack = True, skiprows = 1)
+vis = Re + Im*1j
+
+geom = SourceGeometry(inc= inc, PA= pa, dRA= dra, dDec= ddec)
+
+" Fitting with frank "
+#FF = FrankFitter(r_out, n_pts, geom, alpha = alpha, weights_smooth = w_smooth)
+FF = FourierBesselFitter(r_out, n_pts, geom)
+sol = FF.fit(u, v, vis, Weights)
+#setattr(sol, 'positive', sol.solve_non_negative())
+
+fig = plt.figure(num = 1, figsize = (10, 5))
+plt.plot(sol.r, sol.mean / 1e10)
+plt.xlabel('Radius ["]', size = 15)
+plt.ylabel(r'Brightness profile [$10^{10}$ Jy sr$^{-1}$]', size = 15)
+plt.title('Frank Fit AS209 1mm', size = 15)
+plt.xlim(0, 1.3)
+#plt.savefig('FrankFit_AS209_1mm.jpg', bbox_inches='tight')
+plt.show()
+plt.close()

frank/fourier2d.py

@@ -0,0 +1,100 @@
+
+import numpy as np
+import time
+
+class DiscreteFourierTransform2D(object):
+    def __init__(self, Rmax, N, nu=0):
+        # Remember that now N is to create N**2 points in image plane.
+        self.Xmax = Rmax #Rad
+        self.Ymax = Rmax
+        self.Nx = N 
+        self.Ny = N
+        self.N = self.Nx*self.Ny # Number of points we want to use in the 2D-DFT.
+
+        # Real space collocation points
+        x1n = np.linspace(-self.Xmax, self.Xmax, self.Nx, endpoint=False) # rad
+        x2n = np.linspace(-self.Ymax, self.Ymax, self.Ny, endpoint=False) # rad
+        x1n, x2n = np.meshgrid(x1n, x2n, indexing='ij')
+        # x1n.shape = N**2 X 1, so now, we have N**2 collocation points in the image plane.
+        x1n, x2n = x1n.reshape(-1), x2n.reshape(-1) # x1n = x_array and x2n = y_array
+        self.dx = 2*self.Xmax/self.Nx
+        self.dy = 2*self.Ymax/self.Ny
+
+        # Frequency space collocation points.
+        q1n = np.fft.fftfreq(self.Nx, d = self.dx)
+        q2n = np.fft.fftfreq(self.Ny, d = self.dy)
+        q1n, q2n = np.meshgrid(q1n, q2n, indexing='ij') 
+        # q1n.shape = N**2 X 1, so now, we have N**2 collocation points.
+        q1n, q2n = q1n.reshape(-1), q2n.reshape(-1) # q1n = u_array and q2n = v_array
+
+        self.Xn = x1n
+        self.Yn = x2n
+        self.Un = q1n
+        self.Vn = q2n
+
+    def get_collocation_points(self):        
+        return np.array([self.Xn, self.Yn]), np.array([self.Un, self.Vn])
+
+    def coefficients(self, u = None, v = None, x = None, y = None, direction="forward"):
+        #start_time = time.time()
+        if direction == 'forward':
+            ## Normalization is dx*dy since we the DFT to be an approximation
+            ## of the integral (which depends on the area)
+            norm = 4*self.Xmax*self.Ymax/self.N
+            factor = -2j*np.pi
+
+            X, Y = self.Xn, self.Yn
+            if u is None:
+                u = self.Un
+                v = self.Vn
+        elif direction == 'backward':
+            ## Correcting for the normalization above 1/N is replaced by this:
+            norm = 1 / (4*self.Xmax*self.Ymax)
+            factor = 2j*np.pi
+
+            X, Y = self.Un, self.Vn
+            if u is None:
+                u = self.Xn
+                v = self.Yn
+        else:
+            raise AttributeError("direction must be one of {}"
+                                 "".format(['forward', 'backward']))
+        H = norm * np.exp(factor*(np.outer(u, X) + np.outer(v, Y)))
+        #print("--- %s minutes to calculate 2D-DFT coefficients---" % (time.time()/60 - start_time/60))
+        return H
+
+
+    def transform(self, f, u=None, v=None, direction="forward"):
+        Y = self.coefficients(u, v, direction=direction)
+        return np.dot(Y, f)
+
+
+    @property
+    def size(self):
+        """Number of points used in the 2D-DFT"""
+        return self.N
+
+    @property
+    def uv_points(self):
+        """u and v  collocation points"""
+        return self.Un, self.Vn
+
+    @property
+    def q(self):
+        """Frequency points"""
+        return np.hypot(self.Un, self.Vn)
+
+    @property
+    def Rmax(self):
+        """ Maximum value of the x coordinate in rad"""
+        return self.Xmax
+
+    @property
+    def resolution(self):
+        """ Resolution of the grid in the x coordinate in rad"""
+        return self.dx
+
+    @property
+    def xy_points(self):
+        """ Collocation points in the image plane"""
+        return self.Xn, self.Yn

frank/radial_fitters.py

@@ -28,6 +28,7 @@
 from frank.constants import rad_to_arcsec
 from frank.filter import CriticalFilter
 from frank.hankel import DiscreteHankelTransform
+from frank.fourier2d import DiscreteFourierTransform2D
 from frank.statistical_models import (
     GaussianModel, LogNormalMAPModel, VisibilityMapping
 )
@@ -167,7 +168,7 @@ def interpolate_brightness(self, Rpts, I=None):
         -----
         The resulting brightness will be consistent with higher-resolution fits
         as long as the original fit has sufficient resolution. By sufficient
-        resolution we simply mean that the missing terms in the Fourier-Bessel
+         we simply mean that the missing terms in the Fourier-Bessel
         series are negligible, which will typically be the case if the
         brightness was obtained from a frank fit with 100 points or more.
         """
@@ -429,20 +430,21 @@ class FourierBesselFitter(object):
         R and H should be in arcsec. Assumes a Gaussian vertical structure. 
         Only works with assume_optically_thick=False
     block_size : int, default = 10**5
-        Size of the matrices if blocking is used
+        Size of the matrices if blo cking is used
     verbose : bool, default = False
         Whether to print notification messages
     """
 
-    def __init__(self, Rmax, N, geometry, nu=0, block_data=True,
+    def __init__(self, Rmax, N, geometry=None, nu=0, block_data=True,
                  assume_optically_thick=True, scale_height=None,
-                 block_size=10 ** 5, verbose=True):
+                 block_size=10 ** 5, verbose=True, geometry_on = True):
 
         Rmax /= rad_to_arcsec
 
         self._geometry = geometry
 
         self._DHT = DiscreteHankelTransform(Rmax, N, nu)
+        self._DFT = DiscreteFourierTransform2D(Rmax, N)
 
         if assume_optically_thick:
             if scale_height is not None:
@@ -455,9 +457,10 @@ def __init__(self, Rmax, N, geometry, nu=0, block_data=True,
             model = 'opt_thin'
 
         self._vis_map = VisibilityMapping(self._DHT, geometry, 
-                                          model, scale_height=scale_height,
+                                          model, geometry_on = geometry_on ,scale_height=scale_height,
                                           block_data=block_data, block_size=block_size,
-                                          check_qbounds=False, verbose=verbose)
+                                          check_qbounds=False, verbose=verbose,
+                                          DFT = self._DFT)
 
         self._info  = {'Rmax' : self._DHT.Rmax * rad_to_arcsec,
                        'N' : self._DHT.size
@@ -510,6 +513,8 @@ def _build_matrices(self, mapping):
 
         self._M = mapping['M']
         self._j = mapping['j']
+        self._V = mapping['V']
+        self._Wvalues = mapping['W']
 
         self._H0 = mapping['null_likelihood']
 
@@ -564,7 +569,7 @@ def fit(self, u, v, V, weights=1):
             logging.info('  Fitting for brightness profile using'
                          ' {}'.format(self.fit_method()))
 
-        self._geometry.fit(u, v, V, weights)
+        #self._geometry.fit(u, v, V, weights)
 
         mapping = self.preprocess_visibilities(u, v, V, weights)
         self._build_matrices(mapping)
@@ -574,10 +579,11 @@ def fit(self, u, v, V, weights=1):
     def _fit(self):
         """Fit step. Computes the best fit given the pre-processed data"""
         fit = GaussianModel(self._DHT, self._M, self._j,
-                            noise_likelihood=self._H0)
+                            noise_likelihood=self._H0,
+                            Wvalues= self._Wvalues, V = self._V, DFT = self._DFT)
 
         self._sol = FrankGaussianFit(self._vis_map, fit, self._info,
-                                     geometry=self._geometry.clone())
+                                     geometry=None)
 
         return self._sol