From ff82a845fb185eef75c35650f0cbd4e026aa6b12 Mon Sep 17 00:00:00 2001 From: mariajmellado Date: Wed, 29 May 2024 16:51:58 -0400 Subject: [PATCH] feat: incorporate 2D-DFT to obtain 2D-Frank image --- .../AS209_1mm_FrankFit.py | 0 frank/Examples/DFT_2D_test.ipynb | 195 ++++++++++++++++++ frank/fourier2d.py | 76 +++++++ frank/radial_fitters.py | 7 +- frank/statistical_models.py | 100 ++++----- 5 files changed, 328 insertions(+), 50 deletions(-) rename frank/{fitExample => Examples}/AS209_1mm_FrankFit.py (100%) create mode 100644 frank/Examples/DFT_2D_test.ipynb create mode 100644 frank/fourier2d.py diff --git a/frank/fitExample/AS209_1mm_FrankFit.py b/frank/Examples/AS209_1mm_FrankFit.py similarity index 100% rename from frank/fitExample/AS209_1mm_FrankFit.py rename to frank/Examples/AS209_1mm_FrankFit.py diff --git a/frank/Examples/DFT_2D_test.ipynb b/frank/Examples/DFT_2D_test.ipynb new file mode 100644 index 00000000..7115db73 --- /dev/null +++ b/frank/Examples/DFT_2D_test.ipynb @@ -0,0 +1,195 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "id": "d2808daf-7d75-4451-a75f-d6241049b447", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np \n", + "import matplotlib.pyplot as plt\n", + "\n", + "from frank.geometry import SourceGeometry\n", + "from frank.io import load_uvtable\n", + "from frank.radial_fitters import FrankFitter, FourierBesselFitter" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "ad28005e-bbd1-467f-918d-ba165010fdcd", + "metadata": {}, + "outputs": [], + "source": [ + "# Data source AS209 1mm visibilities gridded\n", + "\n", + "rad_to_arcsec = 3600 * 180 / np.pi\n", + "\n", + "# Huang 2018 \n", + "inc = 34.97\n", + "pa = 85.76\n", + "dra = 1.9e-3\n", + "ddec = -2.5e-3\n", + "r_out = 1.9\n", + "\n", + "# Frank Parameters\n", + "N = 25\n", + "\n", + "# UVtable AS209 at 1mm.\n", + "dir = \"../data/\"\n", + "data_file = dir +'AS209_continuum_prom_1chan_30s_keepflagsFalse_gridded.txt'\n", + "\n", + " # load data\n", + "u, v, Vis, Weights = np.loadtxt(data_file, unpack = True)\n", + "\n", + "geom = SourceGeometry(inc= inc, PA= pa, dRA= dra, dDec= ddec)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "73f358c0-3cf8-4d59-91f3-0521e4531332", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "--- 0.6678717881441116 minutes ---\n" + ] + } + ], + "source": [ + "#FF = FrankFitter(1.9, 300, geom, alpha = alpha, weights_smooth = w_smooth)\n", + "import time\n", + "start_time = time.time()\n", + "FF = FourierBesselFitter(r_out, N, geom)\n", + "sol_new = FF.fit(u, v, Vis, Weights)\n", + "print(\"--- %s minutes ---\" % (time.time()/60 - start_time/60))" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "9282d1fc-9737-4220-a98a-a9d8c4b5527f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "625" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(sol_new.mean) # we corroborate the desire output has 25x25 points." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "4477675b-c2be-4b9f-98e1-f304651ae07d", + "metadata": {}, + "outputs": [], + "source": [ + "dx = dy = 2*r_out/(N**2*rad_to_arcsec)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "f77e0303-e08f-4721-bcb5-00176c97dc13", + "metadata": {}, + "outputs": [], + "source": [ + "I = sol_new.mean.real/(dx*dy)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "c02b56ae-29b5-43fe-8bdd-e7040cafec5f", + "metadata": {}, + "outputs": [], + "source": [ + "def reshape(array, N):\n", + " matrix = np.zeros((N, N), dtype=array.dtype)\n", + "\n", + " # Fill the matrix from bottom to top by columns\n", + " for i in range(N):\n", + " matrix[:, i] = array[i * N:(i + 1) * N][::-1]\n", + " return matrix" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "a6f7e7c6-5819-44e9-a972-1c03be14e169", + "metadata": {}, + "outputs": [], + "source": [ + "I_reshape = reshape(I, N)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "752cdafd-5775-41f6-a7bd-ad6825cdd767", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb8AAAGXCAYAAAAj7F4pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1ZElEQVR4nO3dfXRU9Z0/8PedyczkcQbCQ0JKgCgCChIRBaOCKFkg7eFUsRatZwXququbuGBqPeJWIF1sau22WkXo9lTQPaVQ24Kre5baRh6OFijipuID/oAFiYWEJ/NMZpK5398fMaMjiPeTO3Pnm9z3i3OPZvK5937vPH3y/d7v/VxDKaVARETkIp5UN4CIiMhpTH5EROQ6TH5EROQ6TH5EROQ6TH5EROQ6TH5EROQ6TH5EROQ6TH5EROQ6TH5EROQ6TH5EROQ6TH5ERAmwY8cOzJ07FwUFBTAMA5s3bxat39HRgYULF+Lyyy9HWloabr755nNijh8/jm9961sYM2YMPB4PlixZkpC2uxGTHxFRArS1taG4uBirVq3q1frRaBQZGRn4l3/5F5SWlp43JhwOY8iQIfje976H4uJiO811vbRUN4CIqD8oKytDWVnZF/4+HA7jX//1X/HrX/8ajY2NmDBhAh5//HHMmDEDAJCVlYXVq1cDAN544w00Njaes41Ro0bhqaeeAgA899xzCT8GN2HPj4jIARUVFdi5cyc2bNiAt99+G7fddhvmzJmDAwcOpLpprsTkR0SUZEePHsXatWvx4osvYtq0abj44ovx4IMP4vrrr8fatWtT3TxX4rAnEVGS7du3D9FoFGPGjIl7PBwOY9CgQSlqlbsx+RERJVlrayu8Xi/27t0Lr9cb97vs7OwUtcrdmPyIiJJs0qRJiEajOHHiBKZNm5bq5hCY/IiIEqK1tRUHDx6M/Xz48GHU1tYiNzcXY8aMwZ133om77roL//7v/45Jkybh5MmTqKmpwcSJE/G1r30NAPDee+8hEongzJkzaGlpQW1tLQDgiiuuiG2357HW1lacPHkStbW18Pv9uOyyy5w61H7BUEqpVDeCiKiv27ZtG2688cZzHl+wYAHWrVuHzs5OrFy5Ei+88AL+9re/YfDgwbjmmmtQVVWFyy+/HED3pQwffvjhOdv47Ne0YRjn/H7kyJE4cuRI4g7GBZj8iIjIdXipAxERuQ6THxERuQ4nvBARpVhHRwcikYjt7fj9fqSnpyegRf0fkx8RUQp1dHSgqOgrqK8/Y3tb+fn5OHz4MBOgBUx+REQpFIlEUF9/BkcOb0AwmNnr7TQ3t2NU0e2IRCJMfhYw+RERaSAYzEQwmJXqZrgGkx8RkQ5Ms3uxsz5ZxuRHRKQDJj9H8VIHIiJyHfb8iIh0oFT3Ymd9sozJj4hIB6ayOezJ5CfBYU8iInId9vyIiHTACS+OYvIjItIBk5+jOOxJRESuw54fEZEO2PNzFJMfEZEOlM3kp5j8JDjsSURErtNnkt+qVaswatQopKenY+rUqfjLX/6S6iYlxYoVK2AYRtwybty4VDcrIXbs2IG5c+eioKAAhmFg8+bNcb9XSmHZsmUYNmwYMjIyUFpaigMHDqSmsQnwZce7cOHCc17rOXPmpKaxNlVXV+Pqq69GTk4Ohg4diptvvhkffPBBXExHRwfKy8sxaNAgZGdn49Zbb0VDQ0OKWqwfQ5m2F7KuTyS/jRs3orKyEsuXL8dbb72F4uJizJ49GydOnEh105Ji/PjxOH78eGx5/fXXU92khGhra0NxcTFWrVp13t//6Ec/ws9+9jOsWbMGu3fvRlZWFmbPno2Ojg6HW5oYX3a8ADBnzpy41/rXv/61gy1MnO3bt6O8vBy7du3CH//4R3R2dmLWrFloa2uLxTzwwAN4+eWX8eKLL2L79u04duwY5s2bl8JWa6bnnJ+dhSwzlNK/Js7UqVNx9dVX45lnngEAmKaJwsJC3H///Xj44YdT3LrEWrFiBTZv3oza2tpUNyWpDMPApk2bcPPNNwPo7vUVFBTgO9/5Dh588EEAQFNTE/Ly8rBu3TrcfvvtKWytfZ8/XqC759fY2HhOj7A/OHnyJIYOHYrt27dj+vTpaGpqwpAhQ7B+/Xp84xvfAADs378fl156KXbu3IlrrrkmxS1OnebmZoRCIXz8/9YimGPjfn4t7Rg4ZhGampoQDAYT2ML+SfueXyQSwd69e1FaWhp7zOPxoLS0FDt37kxhy5LnwIEDKCgowEUXXYQ777wTR48eTXWTku7w4cOor6+Pe51DoRCmTp3ab19nANi2bRuGDh2KsWPH4r777sPp06dT3aSEaGpqAgDk5uYCAPbu3YvOzs6413fcuHEYMWJEv359SV/aJ79Tp04hGo0iLy8v7vG8vDzU19enqFXJM3XqVKxbtw5btmzB6tWrcfjwYUybNg0tLS2pblpS9byWbnmdge4hzxdeeAE1NTV4/PHHsX37dpSVlSEajaa6abaYpoklS5bguuuuw4QJEwB0v75+vx8DBgyIi+3Pr68Yhz0dxUsdNFNWVhb7/4kTJ2Lq1KkYOXIkfvOb3+Duu+9OYcso0T47lHv55Zdj4sSJuPjii7Ft2zbMnDkzhS2zp7y8HO+8806/OVftGF7n5yjte36DBw+G1+s9Z1ZYQ0MD8vPzU9Qq5wwYMABjxozBwYMHU92UpOp5Ld36OgPARRddhMGDB/fp17qiogKvvPIKtm7diuHDh8cez8/PRyQSQWNjY1y8m15f0ov2yc/v92Py5MmoqamJPWaaJmpqalBSUpLCljmjtbUVhw4dwrBhw1LdlKQqKipCfn5+3Ovc3NyM3bt3u+J1BoCPPvoIp0+f7pOvtVIKFRUV2LRpE1577TUUFRXF/X7y5Mnw+Xxxr+8HH3yAo0ePuub1/VJKdV+o3utF+7mLWukTw56VlZVYsGABrrrqKkyZMgVPPvkk2trasGjRolQ3LeEefPBBzJ07FyNHjsSxY8ewfPlyeL1e3HHHHalumm2tra1xvZrDhw+jtrYWubm5GDFiBJYsWYKVK1fikksuQVFRER599FEUFBTEzZDsSy50vLm5uaiqqsKtt96K/Px8HDp0CA899BBGjx6N2bNnp7DVvVNeXo7169fjpZdeQk5OTuw8XigUQkZGBkKhEO6++25UVlYiNzcXwWAQ999/P0pKSlw90zMOhz2dpfqIp59+Wo0YMUL5/X41ZcoUtWvXrlQ3KSnmz5+vhg0bpvx+v/rKV76i5s+frw4ePJjqZiXE1q1bFYBzlgULFiillDJNUz366KMqLy9PBQIBNXPmTPXBBx+kttE2XOh429vb1axZs9SQIUOUz+dTI0eOVPfcc4+qr69PdbN75XzHCUCtXbs2FnP27Fn1z//8z2rgwIEqMzNT3XLLLer48eOpa7QmmpqaFAD18b41Knrk+V4vH+9bowCopqamVB9Sn9AnrvMjIuqvYtf5/XU1gjkZvd9Oy1kMLL6P1/lZ1CeGPYmI+j0OezpK+wkvREREicaeHxGRDnhLI0cx+RERacAwTRg2kp+ddd2Iw55EROQ67PkREelAKXsXqnPivkif6fmFw2GsWLEC4XA41U1xhJuO103HCrjreN10rLaxsLWj+sx1fj3XwrjlGhY3Ha+bjhVw1/G66Vh7q+c5atz5YwSzbVzn13oWA0oe5HNtUZ/p+RERESUKz/kREenAVN2LnfXJMu2Sn2maOHbsGHJycmAYRuzx5ubmuP/2d246XjcdK+Cu4+2vx6qUQktLCwoKCuDxJGgAjRVeHKVd8jt27BgKCwu/8PcX+l1/5KbjddOxAu463v56rHV1dXH3LaS+Q7vkl5OTAwD4v8XzkBPwWVqn9ZDsL69TZ7JF8SfaZCehT0f8oviGDusvQ3vU+PKgzzglnGTX0SUbOmnvkm0/KhyZiQqHcrwe2fOTKfwEeGWbR45ftkLI2lv+M/Gy52dowPoLNiyjXbTt4UOaRPGhkZ2ieN/4XFE8Liv68pjPMMePtxzb3NyOolHfjH1fJYSpbPb8OOwpkbTkt2rVKjzxxBOor69HcXExnn76aUyZMuVL1+sZ6swJ+BAMWEsihk+W/DrSZMmpLS0gij8blW0/w2v9G8+E7Ms0IByRMT2yD1CXMF6a/LogWyFNmPz8wnhp8gsIt5/ulW0/wyv7ssz0Wt9BdlpUtO0cn+x9HwzInhtfhuxzCOHMSTOYJds+EHdqxrYUX+f3wx/+EEuXLsXixYvx5JNP2tpWX5CU2Z4bN25EZWUlli9fjrfeegvFxcWYPXs2Tpw4kYzdERGRDXv27MHPf/5zTJw4MdVNcUxSkt9PfvIT3HPPPVi0aBEuu+wyrFmzBpmZmXjuueeSsTsior4vRRe5t7a24s4778QvfvELDBw4MMEHpa+EJ79IJIK9e/eitLT00514PCgtLcXOnTvPiQ+Hw2hubo5biIhcR6lPL3fozdLLYc/y8nJ87Wtfi/vOdoOEn/M7deoUotEo8vLy4h7Py8vD/v37z4mvrq5GVVVVoptBRORKn+9ABAIBBALnP1+6YcMGvPXWW9izZ48TTdNKyiu8LF26FE1NTbGlrq4u1U0iInJegoY9CwsLEQqFYkt1dfV5d1dXV4fFixfjV7/6FdLT0508Ui0kvOc3ePBgeL1eNDQ0xD3e0NCA/Pz8c+Iv9FcJEZFrJOgi97q6urjanl/0/bp3716cOHECV155ZeyxaDSKHTt24JlnnkE4HIZXMDu4r0l4z8/v92Py5MmoqamJPWaaJmpqalBSUpLo3RER9Q92zvd9pjRaMBiMW74o+c2cORP79u1DbW1tbLnqqqtw5513ora2tl8nPiBJ1/lVVlZiwYIFuOqqqzBlyhQ8+eSTaGtrw6JFi5KxOyIiEsrJycGECRPiHsvKysKgQYPOebw/Skrymz9/Pk6ePIlly5ahvr4eV1xxBbZs2XLOJBgiIvqEMrsXO+uTZUmr8FJRUYGKioperx9pUIj4rU3dbWuVnTNsCcviT4VllStOhmVP64mw9SoRLbKKUGjtlE1/buuUfYDao7IqIFHhdGxTGB8VVoQJRGRDO35hEePmTlkFkLYuWXva/bL2dClJ/bRM0bY9J4WvbVR2WdMg/xlRvLTajyfD+qQPT+tZ0bYt0eCuDtu2bbO9jb4i5bM9iYiInKZdYWsiIlfiLY0cxeRHRKQDDYY93YTDnkRE5Drs+RER6YD383MUkx8RkQ447OkoDnsSEZHrsOdHRKQFmxe5g7M9JZj8iIh0wGFPR3HYk4iIXIc9PyIiHbDn5yhtk1/nWQ86u6x1TJvbZbU6PxbW6mwR1ls8E5F1qNsE9To/DsvG9aW1OluisuKhHUoWH4YsvtOIiOKlfKbsvZAhjPdHZR+xTuH2Oyx+RnpElfV6l1FRHVDAMGS1QL0e2Ze1/1CjKD4Y+FgUnxb60HKs0R4WbdsSVnhxlLbJj4jIVdjzcxTP+RERkeuw50dEpAP2/BzF5EdEpAOe83MUhz2JiMh12PMjItKBUt2LnfXJMiY/IiId8JyfozjsSURErsOeHxGRDtjzcxSTHxGRDpTN2Z627gjhPhz2JCIi19G259fR5oPPZ622YFunrAZhY6fssFuE9RNbu0ThOCOo19kcicq2HT0ripfW3mw3WmXbN2Tt6YBs+z6ki+INQ/baeoR/L4bMQaL4SETWfqUyRPEew3qdWq9hvQ4oAPgM2ecwTVoL9ISsZ5OW3iiKz86utxxrdiSh5iyHPR2lbfIjInIVEzaTX8Ja4goc9iQiItdhz4+ISAcc9nQUkx8RkQaUqaBsJDA767oRkx8RkQ5Y3sxRPOdHRESuw54fEZEOeM7PUUx+REQ6YPJzFIc9iYjIddjzIyLSAXt+jmLyIyLSAZOfo7RNfi3t6VBpfkuxbV2ymoIdUdlob0unrMZhu7C2Zzhq/U3bFpXV3pTW6vzYc1IU34mwKL5DNYvipdrN00ndvs8jq0cZ9cjeDFElqwXq7ZS9lw3Deu1QUwnrngrrpKZ7ZZ/bDK+sjmnwlOy96T/aZjk2EpZ9rkg/2iY/IiI3UcrmRe68zk+EyY+ISAcc9nQUZ3sSEZHrsOdHRKQD9vwcxeRHRKQDJj9HcdiTiIhchz0/IiId8K4OjmLyIyLSgDK7Fzvrk3VMfkREOuA5P0fxnB8REbkOkx8RkQ56en52FoHVq1dj4sSJCAaDCAaDKCkpwf/8z/8k6eD0o+2wZyTqQcTwWor9OCI7jFZhbc+IcCz9bJfsTdgejVqODUNWK7LdaJXFq0ZRfIfZJIqPmNbrJ/aGIfx7rss8K4rv9Mjio15ZDchOT4coPgxZfDQy2HKsB7I6pic7ZM99VposPttird9YfHO2KD7zWMRybHtn4ocYnT7nN3z4cPzwhz/EJZdcAqUUnn/+eXz961/H//7v/2L8+PG9b0gfkfCe34oVK2AYRtwybty4RO+GiIhsmDt3Lr761a/ikksuwZgxY/DYY48hOzsbu3btSnXTHJGUnt/48ePxpz/96dOdpGnbwSQi0oOyOeHlk0sdmpvj75wSCAQQCAQuuGo0GsWLL76ItrY2lJSU9L4NfUhSslJaWhry8/OTsWkiov7J/GSxsz6AwsLCuIeXL1+OFStWnHeVffv2oaSkBB0dHcjOzsamTZtw2WWX2WhE35GU5HfgwAEUFBQgPT0dJSUlqK6uxogRI5KxKyIi+oy6ujoEg8HYzxfq9Y0dOxa1tbVoamrCb3/7WyxYsADbt293RQJMePKbOnUq1q1bh7Fjx+L48eOoqqrCtGnT8M477yAnJ+ec+HA4jHD405tOfr7LTkTkBsq0eT+/T9btmb1phd/vx+jRowEAkydPxp49e/DUU0/h5z//ea/b0VckPPmVlZXF/n/ixImYOnUqRo4cid/85je4++67z4mvrq5GVVVVoptBRNS3JGjY01YTTDOuM9KfJf06vwEDBmDMmDE4ePDgeX+/dOlSNDU1xZa6urpkN4mIyPWWLl2KHTt24MiRI9i3bx+WLl2Kbdu24c4770x10xyR9GmYra2tOHToEP7+7//+vL+3MhOJiKjfU58sdtYXOHHiBO666y4cP34coVAIEydOxB/+8Af83d/9nY1G9B0JT34PPvgg5s6di5EjR+LYsWNYvnw5vF4v7rjjjkTvioio30jUOT+rfvnLX/Z6X/1BwpPfRx99hDvuuAOnT5/GkCFDcP3112PXrl0YMmRIondFRETUKwlPfhs2bEjIdjo60+BV1prXacpOXZ61Xk2sl/GyM88tUetllVoNWXmwFuOMKL4jKitX1qVkJ8e7osLyYKZs+1HT+nMJAB5D9hFQXtlrKy2fFk2THa+0PT5YP8WQ3ukTbdvvlZUfOyEshzbAZ63cYY9THbLTKbnNGZZjW7pkbbFEgwkvbsLSK0REGuD9/JzF5EdEpAP2/BzFWxoREZHrsOdHRKQBDns6i8mPiEgHCvaGLhN/i8F+jcOeRETkOuz5ERFpQKnYLfl6vT5Zx+RHRKQBnvNzFoc9iYjIddjzIyLSAa/zcxSTHxGRBjjs6Sxtk19Lpw+mslZbsEsZom2Ho7L49i7ZmeRINHlnnsPoEMVLa2+aqlMUH+lqEcVLa3t2mbLjlVLCbwxp7VCvR1bv0hCeiZDWJm33WLvDNwA0m5mibfvCsrZnpclqhzYIa4FmeWXPzak267U927p4xqiv0zb5ERG5CWd7OovJj4hIB6bRvdhZnyxj352IiFyHPT8iIg1wwouzmPyIiDSglAElnLz3+fXJOiY/IiINsOfnLJ7zIyIi12HPj4hIA0rZ7PnxUgcRJj8iIg3wnJ+zOOxJRESuw54fEZEOTAOKF7k7Rtvk5zG6FyuEpTfl8cJx+E5TtoMwuqxv25DW6oyK4ruUrHZl1JS2x/qxAoAprKWphNuXlsL3eNKTuXl0RWW1TKNe2fPfbjRajm0xrNe6BICAKfs6aemUxZ8V1uRt7JQNbLV0Wq812ib9UrCA5c2cxWFPIiJyHW17fkREbsIJL85i8iMi0oCyec7P1vlCF+KwJxERuQ57fkREGuCEF2cx+RERaYDn/JzFYU8iInId9vyIiDRgmgZMG5NW7KzrlBdeeCFh27rrrrtsrc/kR0SkATec81u4cCEMIzFJmsmPiKgfcMs5v+LiYnz961/v9fqbN2/G22+/bbsdTH5EROSYK664AsuXL+/1+keOHOnfya8t6oWC11Jsh7Dmn7D0JqLieNkKSlAA0hQWi+xSstqPSlgLVAlvQNYVbZdt35TVulSQtd8wAqL4aLRNFK88sufHMGRz0Dq6mkTxaZ5My7FnPbJj7VQDRPER4QexMSKLD/pkz2WjoLZnezQZtT37f88vGAwiM9P6e/B8MjIyEAwGbbdF2+RHROQmpjJg2khgdtZ1SmNjo+1tPPvss3j22Wdtb4eXOhARkeuw50dEpAHW9nQWkx8RkQbccKmDTjjsSUREKbVlyxZcdNFFju6TPT8iIg2YsDnhBX132LOtrQ0ffviho/tk8iMi0kB/vNRh2bJlluLef//9JLfkXEx+REQuVF1djd///vfYv38/MjIycO211+Lxxx/H2LFjE7aPlStXYsCAAQiFQheMa2+XXf+bCEx+REQaUDav85P2/LZv347y8nJcffXV6OrqwiOPPIJZs2bhvffeQ1ZWVq/b8VkXX3wxpk2bhueee+6Ccb/97W8xf/78hOzTKiY/IiINOD3suWXLlrif161bh6FDh2Lv3r2YPn16r9vxWSUlJfjzn//8pXGGYUA5PF2VyY+ISAPmJ4ud9QGgubk57vFAIIBA4MvL+DU1dZfKy83NtdGKeIsXL8Ybb7zxpXE33HADtm7dmrD9WtEvkl+n8C8eaa3OsHAFE8n7C6YLslqdpuoUxneJ4pNOWOtSerGTUrLaobBYb/bT7cueT2mtVNOUbb/LtH5uJeqRvXdaDVkt0LaudFF8TtQvio8IM0lrl/X32tmovleJFRYWxv28fPlyrFix4oLrmKaJJUuW4LrrrsOECRMS1pbJkydj8uTJXxo3ePBg3HDDDQnbrxX9IvkREfV1iRr2rKuriyv8bKXXV15ejnfeeQevv/56r/ff14j/fNmxYwfmzp2LgoICGIaBzZs3x/1eKYVly5Zh2LBhyMjIQGlpKQ4cOJCo9hIR9Uum+rS4de+W7u0Eg8G45cuSX0VFBV555RVs3boVw4cPd+BIuxP0TTfd5Mi+vog4+bW1taG4uBirVq067+9/9KMf4Wc/+xnWrFmD3bt3IysrC7Nnz0ZHh3R4iYiIkkUphYqKCmzatAmvvfYaioqKHNt3e3s7tm/f7tj+zkc87FlWVoaysrLz/k4phSeffBLf+973YnfqfeGFF5CXl4fNmzfj9ttvt9daIqJ+yunZnuXl5Vi/fj1eeukl5OTkoL6+HgAQCoWQkZHR63b0FQk9a3v48GHU19ejtLQ09lgoFMLUqVOxc+fO864TDofR3NwctxARuU33sKe9RWL16tVoamrCjBkzMGzYsNiycePG5BygZhI64aXnL4e8vLy4x/Py8mK/+7zq6mpUVVUlshlERPQlnL6uTjcpn6+7dOlSNDU1xZa6urpUN4mIyHE9w552FrIuoT2//Px8AEBDQwOGDRsWe7yhoQFXXHHFedexegEmEVF/ZsKwdWeGvnxXh1RIaM+vqKgI+fn5qKmpiT3W3NyM3bt3o6SkJJG7IiIi6jVxz6+1tRUHDx6M/Xz48GHU1tYiNzcXI0aMwJIlS7By5UpccsklKCoqwqOPPoqCggLcfPPNiWw3EVG/4rY7uaf6nKM4+b355pu48cYbYz9XVlYCABYsWIB169bhoYceQltbG/7xH/8RjY2NuP7667Flyxakp8tKGRERuUnPxep21u8r8vPzsXr16pS2QZz8ZsyYccGMbRgGvv/97+P73/++rYZJSMdupbU9vUl+T3Uhajk2DbLzo4awNqa0tqTH8AnjZW+5qIqI4uWlgZP74hrC45W+XtJ4U/D8dKqzSds2AESFf/mbwvj2LtlrK9l8X+tl6SYUCuGf/umfUtqGlM/2JCIiQH0y4aW3i+qDE14eeOABrF+/PiX7ZvIjItJAzzk/O0tf88wzz+BXv/pVSvbNuzoQEWnATef8egwfPhw5OTkp2Td7fkRElBK33HILduzYgUhEem7fPiY/IiINqE/O29lZ+ppHH30UXq8XixcvdnzfTH5ERBpwurC1DubNm4fi4mL8x3/8B6ZPn47du3c7tm+e8yMiopT47D39Xn/9dVx77bXIz8/H9OnTMWnSpNgyePDghO+byY+ISANunPDSUyHsr3/9K2pra1FbW4sjR45g48aN2LhxIwyj+5gKCgpw5ZVX4qWXXkrYvpn8iIg0YPe8XV885zdy5EiMHDkydvNzoLsedE8y7Pnvu+++i1deeSWh+2byIyKilHjvvffQ0NCAKVOmICsrCwAQDAYxbdo0TJs2LRYXjUaxf//+hO6bE16IiDTgxgkvTzzxBEpLS/Hee+/FPd7Q0IAf/OAHWLlyJd5++214vV6MHz8+ofvWtucXVQaiFsewpbU6kz044EniHnzC2p5eYe1Naa3IqBkWxUNcu9Iv276UsJapEtRh7SbbvmnKrndSngxRfDJJn5tOJYwXlm2V5oKzpvX35lkz8Z9xNw577ty5E6NHj8bVV18deywcDqOkpAQffvghlFJYsWIFqqur8d3vfjeh+2bPj4iIUuL48eMYM2ZM3GMbNmzAkSNHcNVVV+GnP/0pLr74Yjz88MN44403ErpvJj8iIg24cdgzHA4jOzs77rHf/e538Hq92LhxIxYvXow//elPSEtLw1NPPZXQfWs77ElE5CZuvNThK1/5Co4cORL7ub29HTU1Nbj22msxatQoAEBhYSGmTZvGnh8REfUPM2bMwJ49e/D2228DAF544QWcPXsWZWVlcXH5+fk4depUQvfN5EdEpAGVgKWv+e53vwufz4ebbroJt9xyCyorK+H1ejF//vy4uNOnTyMYDCZ030x+REQaUDBiQ5+9WfribM9x48Zh06ZNSE9Px0svvYRwOIyqqioUFRXFYkzTxJ49ezB8+PCE7pvn/IiINGBCemHMuev3RXPmzMHRo0dx4MABhEIh5Ofnx/3+1VdfxZkzZ87pDdrF5EdERCnl8XgwduzYL/zdokWLMG/evITuk8mPiEgDShlQNmZs2llXZ7NmzcKsWbMSvl0mPyIiDbh12DNVOOGFiIhcR9ueXzIrFng0Gx3wwXr9zSg6Rdv2CLYNAF6PsHaoR1aL0ox2ieKlf88ahuwtrSBrjyfZtUaFtU99aZmieI/g7900Q/heULL3mrQGblTJvhD6WsUTu995fe14U03b5EdE5CZuLGydShz2JCIi12HPj4hIAxz2dBZ7fkREGugZ9rSz9DVVVVX46KOPUrJvJj8iIkqJnlJmc+fOxX/913/BNJ27YIPJj4hIA268n9/KlSsxYsQI/Pd//zduueUWFBYW4tFHH427zVGyMPkREWnAjcnvkUcewaFDh/Dqq6/itttuw+nTp/HYY49h9OjRmDNnDn73u9+hq0t6eZQ1TH5ERJRSpaWl2LBhA/72t7/hxz/+McaOHYtXX30V3/zmNzF8+HA8/PDDOHDgQEL3yeRHRKQBN054+bxBgwahsrIS7777Ll5//XXccccdOHHiBJ544gmMGzcOM2fOxKZNmxKyLyY/IiINKJtDnsICOFo7dOgQXn75ZdTU1MQeGz58OLZu3YpvfOMbmDJlCurq6mztg8mPiEgDZgKWvqyzsxMbNmzAzJkzMWbMGDz++OPo6upCZWUl9u/fjw8//BBvvPEGysrK8Oabb6KiosLW/rS9yF3SjTeEvX3piWFpLdA04QpKML3Xp4T1Fg1ZvcWAN0cUb5qyWqMKUVG8Iax1aSrZyXEljJdK88hqb8prq8ri/d5sy7GG8G9jn7C2pyH84PqTXJQ3KvhekMTShb3//vv4xS9+gf/8z//EmTNnoJTCtddei3vvvRe33XYbAoFP3+MlJSV45ZVXcM0112D79u229qtt8iMichM33s/v+uuvx86dO6GUQjAYxH333Yd7770XEyZMuOB648ePx549e2ztm8mPiEgDbryf35///GdceeWVuPfee/Gtb30LmZnWRkr+4R/+AdOnT7e1byY/IiJKiT179mDy5Mni9UpKSlBSUmJr35zwQkSkAacvct+xYwfmzp2LgoICGIaBzZs3J+W4LqQ3iS9RmPyIiDSgErBItLW1obi4GKtWrUpI+/saDnsSEblQWVkZysrKUt2MlGHyIyLSQPfQZe9nbPYMezY3N8c9HggE4i4XSCWv19ur9QzDQFZWFkaMGIEbb7wRlZWVGDVqlK22cNiTiEgDiRr2LCwsRCgUii3V1dWOHseFKKV6tZimiZaWFrz77rt45plnMGnSJOzbt89WW9jzIyLqR+rq6hAMBmM/69LrA2Drfn3hcBj/93//h6effhpr1qzB8uXL8fvf/77X22PyIyLSgN3bEvWsGwwG45JffxEIBHDppZfi2WefxRtvvIEdO3bY2h6HPYmINOD22p4SEyZMQGNjo61taNvz8xkKPsPan0FW43pIa2/6vbL4gEf2N0UmrA9LhJEu2naaYNsAkGYIt+/NEMUr4UdUqeR+pA3Dn9TtS58fX1qWKD6Zz49P+t6BbDKDT1i3VVra0yf8014S3h96Da2trTh48GDs58OHD6O2tha5ubkYMWJEUvZ55swZpKenW67kcj7t7e14+umn8YMf/MBWW8Sv4ZddGLlw4UIYhhG3zJkzx1YjiYj6O6XsLxJvvvkmJk2ahEmTJgEAKisrMWnSJCxbtiwJR9dtyJAhuP/++21to7y8HEOHDsXIkSNtbUfc8+u5MPLb3/425s2bd96YOXPmYO3atbGfdTrhSkSkIwUDpo0b0kpvZjtjxgwoh28C2DN7MxHbsUuc/KxcGBkIBJCfn9/rRhERuU1vem+fX78veP311/Htb3/b1vqJkJRzftu2bcPQoUMxcOBA3HTTTVi5ciUGDRqUjF0REVEfcvDgwbhzjb0hvRfk+SQ8+c2ZMwfz5s1DUVERDh06hEceeQRlZWXYuXPnea/uD4fDCIfDsZ8/X52AiMgN3HBLo61bt6a6CTEJT36333577P8vv/xyTJw4ERdffDG2bduGmTNnnhNfXV2NqqqqRDeDiKhPSdR1fjq74YYbUt2EmKTP2L3oooswePDgL+zmLl26FE1NTbGlrq4u2U0iIiKXS/p1fh999BFOnz6NYcOGnff3OhVdJSJKld7clujz65N14uR3oQsjc3NzUVVVhVtvvRX5+fk4dOgQHnroIYwePRqzZ89OaMOJiPoTNwx76kSc/N58803ceOONsZ8rKysBAAsWLMDq1avx9ttv4/nnn0djYyMKCgowa9Ys/Nu//Rt7d0REpA1x8vuyCyP/8Ic/2GoQEZEbueU6P11oW9uzS3UvVghLbyIgnObjkxYVFPIZ1msipitZrchMCKu7Cw+1y9MhijdVlyje65eNGCgVFcV3RWXtT/PKap9Ka2/6PMJaoML4DGOg9W0r2XOfDlmd1EzhjU3ThR906feC5NKxBFxmdg43XOqgk/5Qn5WIiEhE254fEZGbcMKLs5j8iIg04IZLHc5X5csqwzDQ1SU7bXIhTH5EROQIO3djSPQdKJj8iIg04IZhT9PUZ1oOkx8RkQYUDPE9+T6/PlnH5EdEpAEFe723PtDx0wovdSAiItdhz4+ISANuOOenEyY/IiINuOFSB51w2JOIiFxH255fljeKrDRrdRobO5N7GD7hJKqMNNkK7VHrf4MEorJjzVTZonjTkE1FVh5h7UojUxTfqdpF8VEzLIr3emT1K9MMWbxhyP6+TDNktUPTDVnt1mw1wHKsX8naEjBk782AV/jcCGvsZgm/FoJp1t/LPiPx/SwOezpL2+RHROQm6pN/dtYn6zjsSURErsOeHxGRBjjs6SwmPyIiDXC2p7M47ElERK7Dnh8RkQY47OksJj8iIg0o1b3YWZ+sY/IjItKA+cliZ32yjuf8iIjIddjzIyLSAM/5OYvJj4hIBzbP+fFaBxltk58JA1FlrZafR/iqB7yyeL9XVlNQWtszvctrOTaqfKJtK+Gfg6YpO3PgEdau9Bqy9gcMWW3SqKdTFN+lZLVAPbD+WgGAR1jv0g9Z7dNMJavtmaEyLMcGYT0WALK8smNNF36uMoXfVj6P7L2f7rX+3jd5hq3P0zb5ERG5CSe8OIvJj4hIA7zUwVmc7UlERK7Dnh8RkQY47OksJj8iIg0opaBsjF3aWdeNOOxJRESuw54fEZEGeJG7s5j8iIg0wPv5OYvDnkRE5Drs+RERaYDDns7SNvkZnyxW+KVljIT93QxZRSu0d8nKNmUJyqEpJWuMKSyHZqgsUbx0rMWHgCg+bJwVxSvhhG/TiIrivZA9n1IBQfkxAMhUsnJoGYL2D/D5ZW1Jk32wsnyyz0m28NsqlCZ7L2R6Je8F2fvGCiY/Z3HYk4hIA93n/Oz8651Vq1Zh1KhRSE9Px9SpU/GXv/wlkYelLSY/IiKX2rhxIyorK7F8+XK89dZbKC4uxuzZs3HixIlUNy3pmPyIiDTQM+xpZ5H6yU9+gnvuuQeLFi3CZZddhjVr1iAzMxPPPfdc4g9QM0x+REQa6ClsbWeRiEQi2Lt3L0pLS2OPeTwelJaWYufOnQk+Ov1oO+GFiIjkmpub434OBAIIBM6daHbq1ClEo1Hk5eXFPZ6Xl4f9+/cntY06YM+PiEgDCgqmjaVnykthYSFCoVBsqa6uTvGR6Yk9PyIiDSTqfn51dXUIBoOxx8/X6wOAwYMHw+v1oqGhIe7xhoYG5Ofn974hfQR7fkRE/UgwGIxbvij5+f1+TJ48GTU1NbHHTNNETU0NSkpKnGpuyrDnR0SkgVTcz6+yshILFizAVVddhSlTpuDJJ59EW1sbFi1aZKMlfQOTHxGRBlJxP7/58+fj5MmTWLZsGerr63HFFVdgy5Yt50yC6Y+Y/IiIXKyiogIVFRWpbobjtE1+aYYJn8daRz7dK+vw+6OyU51eQ1aDMF1YCzSqrG/fVLK2m0r2EntM2bF6zGxRfKeS1UTsVOmi+C5hzUWP8LS3x3LF2W5eyN4MAeFH0m/Itp/utR6f7FqdOcIyqTk+Wc9GWvM3K61LEC2JtYa1PZ2lbfIjInKTnksW7KxP1on+tKuursbVV1+NnJwcDB06FDfffDM++OCDuJiOjg6Ul5dj0KBByM7Oxq233nrOVFoiIqJUEiW/7du3o7y8HLt27cIf//hHdHZ2YtasWWhra4vFPPDAA3j55Zfx4osvYvv27Th27BjmzZuX8IYTEfUnCjbLm6X6APoY0bDnli1b4n5et24dhg4dir1792L69OloamrCL3/5S6xfvx433XQTAGDt2rW49NJLsWvXLlxzzTWJazkRUT/CYU9n2brIvampCQCQm5sLANi7dy86OzvjCqWOGzcOI0aM+MJCqeFwGM3NzXELEZHbOF3Y2u16nfxM08SSJUtw3XXXYcKECQCA+vp6+P1+DBgwIC42Ly8P9fX1591OdXV1XB26wsLC3jaJiIjIkl4nv/LycrzzzjvYsGGDrQYsXboUTU1NsaWurs7W9oiI+iI7Ra3tDpm6Ua8udaioqMArr7yCHTt2YPjw4bHH8/PzEYlE0NjYGNf7u1Ch1C+63QYRkZuYyuY5P457ioh6fkopVFRUYNOmTXjttddQVFQU9/vJkyfD5/PFFUr94IMPcPToUVcUSiUior5B1PMrLy/H+vXr8dJLLyEnJyd2Hi8UCiEjIwOhUAh33303KisrkZubi2AwiPvvvx8lJSWc6UlEdAHqM/fk6+36ZJ0o+a1evRoAMGPGjLjH165di4ULFwIAfvrTn8Lj8eDWW29FOBzG7Nmz8eyzzyaksURE/ZWCvbs6MPXJiJKflarh6enpWLVqFVatWtXrRgGAz2O9tqdPWMMvzZDFp3tl8V2CWp2ArLandFhfCec0eaPCOqamrLZke1RWE1HBL4qPKtnXh/Qci09YS9PvkT3/PkMWn54me72yfda3n+6VbXuA7KVCUFirM1tYw3egX/Zey/F1Wo41DOuxpCfW9iQi0gAvcncWkx8RkQaUsnnOj7M9RWxVeCEiIuqL2PMjItIAhz2dxeRHRKQBJj9ncdiTiIhchz0/IiIN9FTotLM+WcfkR0SkAQ57OovJj4hIA0x+zuI5PyIich32/IiINGB+8s/O+mSdtsnP7zHht1rbM8m1OtOF76mocPTBFMT7PbJ6i9KPg3T70oGWdGHt0KjwAKTHK62KkSZ8fnzCeGE5TXH9zaDfeny6rIwpBvqTW6tzUEBWq3OgPyKKD2WELcd6umTbtkIZCsqwM+GFw54SHPYkIiLX0bbnR0TkJsrmhBf2/GSY/IiINGDChMFzfo7hsCcREbkOe35ERBpghRdnMfkREWnANEwYNmZ7cthThsOeRETkOuz5ERFpgBNenMXkR0SkASY/Z3HYk4iIXIc9PyIiDXC2p7O0TX4eQ8FjsWanz2IN0E/jk1uDUClZh1pS/rFDWBsTkMVHhJ+fTmG815A9N4b0cIWEpT2RJhwrCQhre/qF9TR9wvZkCLYf9MmenEzh52SwX1arc4CwVufAjA5RfCh41nKspzPxtT1NRGEgamt9sk7b5EdE5CbqkwJndtYn63jOj4iIXIc9PyIiDfAid2cx+RERaaD7nF/vB+N4zk+Gw55EROQ67PkREWnB3qUO4LCnCJMfEZEGTBWFncG47vXJKg57EhGR67DnR0SkAVZ4cRZ7fkREGlCI2l6S5bHHHsO1116LzMxMDBgwIGn7cRKTHxERXVAkEsFtt92G++67L9VNSRhthz3TPKblmp3S2p5+YbwprLcYFVcZsv43iNdivdMePmFtyS7hyEmX8Fg7TVl7TOH2pbUupbVJ04XvBamsNNkBB4THK6m/meGVtSXkk9XqzA2ERfHBgKye5oAc67U6ASA9ZL39kYjsWK3ovkhdz4vcq6qqAADr1q1L2j6cpm3yIyJyk0TV9mxubo57PBAIIBAI2Gpbf8RhTyKifqSwsBChUCi2VFdXp7pJWmLPj4hIA0pFoYS3IPv8+gBQV1eHYDAYe/yLen0PP/wwHn/88Qtu8/3338e4ceN63SadMfkREWkgUef8gsFgXPL7It/5znewcOHCC8ZcdNFFvW6P7pj8iIg00H25go2en/BShyFDhmDIkCG93l9fx+RHREQXdPToUZw5cwZHjx5FNBpFbW0tAGD06NHIzs5ObeN6icmPiEgDStms8KKSd6nDsmXL8Pzzz8d+njRpEgBg69atmDFjRtL2m0yc7UlEpAEzAf+SZd26dVBKnbP01cQHMPkREZELcdiTiEgDibrUgazRLvkp1V2loK3LeimjNmFNrvaoLL4jKusgnxXHW3/DS8uJRYTlxJJd3kwaLy1vJh34kZY3SzZp+Trp82PA+gpKEAsAPo+s5FegS1bezOPtFMWndcrKoZkR6/Etke629HxfJUKiKryQNdolv5aWFgDAV3f/MsUtISK6sJaWFoRCoVQ3g3pBu+RXUFCAuro65OTkwDA+7bU0NzejsLDwnOoF/ZWbjtdNxwq463j767EqpdDS0oKCgoIEbtO0Oeyp2TCG5rRLfh6PB8OHD//C31utXtBfuOl43XSsgLuOtz8ea+J7fFGbA5c85yfB2Z5EROQ62vX8iIjcqHvYksOeTukzyS8QCGD58uWuuS+Vm47XTccKuOt43XSsdjH5OctQiZyrS0REIs3NzQiFQsjJHAfD8PZ6O0pF0dK+H01NTf3u/Goy8JwfERG5Tp8Z9iQi6s847OksJj8iIg3YLU/G8mYyHPYkIiLXYc+PiEgD3bU5WdvTKUx+REQasHvOjuf8ZDjsSURErsOeHxGRBtjzcxaTHxGRBuzcyy8R67sNhz2JiMh12PMjItIAhz2dxeRHRKQBJj9ncdiTiIhchz0/IiIt2O25secnweRHRKQBDns6i8mPiEgDvNTBWTznR0RErsOeHxGRBpSyWdhasbC1BJMfEZEWorBzM1vwrg4iHPYkIiLXYc+PiEgD3bM1e9/z47CnDJMfEZEW7CU/DnvKcNiTiIhchz0/IiId2Bz2BIc9RZj8iIg0oGwOW9pd32047ElERK7Dnh8RkRY44cVJTH5ERFpQNvMXk58Ehz2JiMh12PMjItKC4qQVB7HnR0SUQn6/H/n5+eiu7Wlvyc/Ph9/vd/wY+iJDsSYOEVFKdXR0IBKJ2N6O3+9Henp6AlrU/zH5ERGR63DYk4iIXIfJj4iIXIfJj4iIXIfJj4iIXIfJj4iIXIfJj4iIXIfJj4iIXOf/A7yWd2hNcKdIAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(15, 15))\n", + "plt.matshow(I_reshape, cmap=\"magma\")\n", + "cmap = plt.colorbar()\n", + "cmap.set_label(r'I [Jy $sr^{-1}$]', size = 15)\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "PROD_Frank2DVenv", + "language": "python", + "name": "prod_frank2dvenv" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/frank/fourier2d.py b/frank/fourier2d.py new file mode 100644 index 00000000..42b85290 --- /dev/null +++ b/frank/fourier2d.py @@ -0,0 +1,76 @@ + +import numpy as np + +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 = 2*Rmax # [Rmax] = rad. + self.Ymax = 2*Rmax + self.Nx = N + self.Ny= N + self.N = N**2 # Number of points we want to use in the 2D-DFT. + + # Real space collocation points + x1n = np.linspace(0, self.Xmax, self.Nx) # rad + x2n = np.linspace(0, self.Ymax, self.Ny) # rad + x1n, x2n = np.meshgrid(x1n, x1n, 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 + dx = 2*self.Xmax/self.N + dy = 2*self.Ymax/self.N + + # Frequency space collocation points. + # The [1:] is because to not consider the 0 baseline. But we're missing points. + q1n = np.fft.fftfreq(self.Nx+1, d = dx)[1:] + q2n = np.fft.fftfreq(self.Ny+1, d = dy)[1:] + 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, direction="forward"): + if direction == 'forward': + norm = 1 + elif direction == 'backward': + norm = 1 / self.N + else: + raise AttributeError("direction must be one of {}" + "".format(['forward', 'backward'])) + if u is None: + u = self.Un + v = self.Vn + + factor = -2j*np.pi/self.N + if direction == "forward": + H = norm * np.exp(factor*(np.outer(u, self.Xn) + np.outer(v, self.Yn))) + else: + H = norm * np.matrix(np.exp(factor*(np.outer(u, self.Xn) + np.outer(v, self.Yn)))).getH() + return H + + + def transform(self, f, u, v, 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) \ No newline at end of file diff --git a/frank/radial_fitters.py b/frank/radial_fitters.py index dfb6f87f..0f606874 100644 --- a/frank/radial_fitters.py +++ b/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 ) @@ -443,6 +444,7 @@ def __init__(self, Rmax, N, geometry, nu=0, block_data=True, self._geometry = geometry self._DHT = DiscreteHankelTransform(Rmax, N, nu) + self._DFT = DiscreteFourierTransform2D(Rmax, N) if assume_optically_thick: if scale_height is not None: @@ -457,7 +459,8 @@ def __init__(self, Rmax, N, geometry, nu=0, block_data=True, self._vis_map = VisibilityMapping(self._DHT, geometry, model, 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 @@ -577,7 +580,7 @@ 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, - Wvalues= self._Wvalues, V = self._V) + Wvalues= self._Wvalues, V = self._V, DFT = self._DFT) self._sol = FrankGaussianFit(self._vis_map, fit, self._info, geometry=self._geometry.clone()) diff --git a/frank/statistical_models.py b/frank/statistical_models.py index 183cdbcd..a1a61f61 100644 --- a/frank/statistical_models.py +++ b/frank/statistical_models.py @@ -66,7 +66,8 @@ class VisibilityMapping: """ def __init__(self, DHT, geometry, vis_model='opt_thick', scale_height=None, block_data=True, - block_size=10 ** 5, check_qbounds=True, verbose=True): + block_size=10 ** 5, check_qbounds=True, verbose=True, + DFT = None): _vis_models = ['opt_thick', 'opt_thin', 'debris'] if vis_model not in _vis_models: @@ -81,6 +82,7 @@ def __init__(self, DHT, geometry, self._chunk_size = block_size self._DHT = DHT + self._DFT = DFT self._geometry = geometry # Check for consistency and report the model choice. @@ -178,8 +180,8 @@ def map_visibilities(self, u, v, V, weights, frequencies=None, geometry=None): frequencies = np.ones_like(V) channels = np.unique(frequencies) - Ms = np.zeros([len(channels), self.size, self.size], dtype='f8') - js = np.zeros([len(channels), self.size], dtype='f8') + Ms = np.zeros([len(channels), self.size, self.size], dtype='c8') + js = np.zeros([len(channels), self.size], dtype='c8') for i, f in enumerate(channels): idx = frequencies == f @@ -199,11 +201,13 @@ def map_visibilities(self, u, v, V, weights, frequencies=None, geometry=None): Ndata = len(Vi) while start < Ndata: qs = qi[start:end] + us = u[start:end] + vs = v[start:end] ks = ki[start:end] ws = wi[start:end] Vs = Vi[start:end] - X = self._get_mapping_coefficients(qs, ks) + X = self._get_mapping_coefficients(qs, ks, us, vs) wXT = np.array(X.T * ws, order='C') @@ -462,7 +466,7 @@ def interpolate(self, f, r, space='Real'): return self._DHT.interpolate(f, r, space) - def _get_mapping_coefficients(self, qs, ks, geometry=None, inverse=False): + def _get_mapping_coefficients(self, qs, ks, u, v, geometry=None, inverse=False): """Get :math:`H(q)`, such that :math:`V(q) = H(q) I_\nu`""" if self._vis_model == 'opt_thick': @@ -486,7 +490,8 @@ def _get_mapping_coefficients(self, qs, ks, geometry=None, inverse=False): else: direction='forward' - H = self._DHT.coefficients(qs, direction=direction) * scale + #H = self._DHT.coefficients(qs, direction=direction) * scale + H = self._DFT.coefficients(u, v, direction=direction) * scale return H @@ -529,7 +534,8 @@ def Rmax(self): @property def q(self): r"""Frequency points, unit = :math:`\lambda`""" - return self._DHT.q + #return self._DHT.q + return self._DFT.q @property def Qmax(self): @@ -539,7 +545,8 @@ def Qmax(self): @property def size(self): """Number of points in reconstruction""" - return self._DHT.size + #return self._DHT.size + return self._DFT.size @property def scale_height(self): @@ -631,9 +638,10 @@ class GaussianModel: def __init__(self, DHT, M, j, p=None, scale=None, guess=None, Nfields=None, noise_likelihood=0, - Wvalues = None, V = None): + Wvalues = None, V = None, DFT = None): self._DHT = DHT + self._DFT = DFT self._Wvalues = Wvalues self._V = V @@ -691,24 +699,13 @@ def __init__(self, DHT, M, j, p=None, scale=None, guess=None, en = (n+1)*Nr self._Sinv[sn:en, sn:en] += Sj[n] - else: + else: # p is None, so we are interested in this case. + " New GP Kernel below" #self._Sinv = None - q_array = self._DHT.q - - def true_squared_exponential_kernel(q, p, l): - - q1, q2 = np.meshgrid(q, q) - p1, p2 = np.meshgrid(p, p) - SE_Kernel = np.sqrt(p1 * p2) * np.exp(-0.5*(q1-q2)**2 / l**2) - return SE_Kernel - - Ykm = self._DHT.coefficients(direction="backward") - # We continue after set M matrix because is needed to calculate - # best parameters for S matrix. # Compute the design matrix - self._M = np.zeros([Nr*Nfields, Nr*Nfields], dtype='f8') - self._j = np.zeros(Nr*Nfields, dtype='f8') + self._M = np.zeros([Nr*Nfields, Nr*Nfields], dtype='c8') + self._j = np.zeros(Nr*Nfields, dtype='c8') for si, Mi, ji in zip(self._scale, M, j): for n in range(0, Nfields): @@ -724,35 +721,42 @@ def true_squared_exponential_kernel(q, p, l): self._like_noise = noise_likelihood - # M is already defined, so we find best parameters for S matrix and use it. - m, c , l = self.minimizeS() - pI = np.exp(m*np.log(q_array) + c) - S_fspace = true_squared_exponential_kernel(q_array, pI, l) - S_real = np.dot(np.transpose(Ykm), np.dot(S_fspace, Ykm)) + " New GP " + self.u, self.v = self._DFT.uv_points + self.Ykm = self._DFT.coefficients(direction="backward") + self.Ykm_f = self._DFT.coefficients(direction="forward") + m, c , l = -4.8, 59.21, 1.21e5 + #m, c, l = self.minimizeS() # Finding best parameters to S matrix. + S_real = self.calculate_S(self.u, self.v, l, m, c) S_real_inv = np.linalg.inv(S_real) self._Sinv = S_real_inv self._fit() + def true_squared_exponential_kernel(self, u, v, l, m, c): + u1, u2 = np.meshgrid(u, u) + v1, v2 = np.meshgrid(v, v) + q1 = np.hypot(u1, v1) + q2 = np.hypot(u2, v2) + + def power_spectrum(q, m, c): + return np.exp(m*np.log(q) + c) + p1 = power_spectrum(q1, m, c) + p2 = power_spectrum(q2, m, c) + + SE_Kernel = np.sqrt(p1 * p2) * np.exp(-0.5*((u1-u2)**2 + (v1-v2)**2)/ l**2) + return SE_Kernel + + def calculate_S(self, u, v, l, m, c): + S_fspace = self.true_squared_exponential_kernel(u, v, l, m, c) + S_real = np.matmul(self.Ykm, np.matmul(S_fspace, self.Ykm_f)) + return S_real + def minimizeS(self): from scipy.optimize import minimize from scipy.special import gamma V = self._V - def calculate_S(m, c, l): - q_array = self._DHT.q - p_array = c*(q_array**m) - def true_squared_exponential_kernel(q, p, l): - q1, q2 = np.meshgrid(q, q) - p1, p2 = np.meshgrid(p, p) - SE_Kernel = np.sqrt(p1 * p2) * np.exp(-0.5*(q1-q2)**2 / l**2) - return SE_Kernel - - Ykm = self._DHT.coefficients(direction="backward") - S_fspace = true_squared_exponential_kernel(q_array, p_array, l) - S_real = np.dot(np.transpose(Ykm), np.dot(S_fspace, Ykm)) - return S_real - def calculate_D(S): S_real_inv = np.linalg.inv(S) Dinv = self._M + S_real_inv @@ -782,7 +786,7 @@ def likelihood(param, data): def inv_gamma_function(l, alpha, beta): return ((gamma(alpha)*beta)**(-1))*((beta/l)**(alpha + 1))*np.exp(-beta/l) - S = calculate_S(m,c, l) + S = self.calculate_S(self.u, self.v, l, m, c) [Dinv, D] = calculate_D(S) mu = calculate_mu(Dinv) logdetS = np.linalg.slogdet(S)[1] @@ -796,7 +800,7 @@ def inv_gamma_function(l, alpha, beta): - 0.5*(factor + logdetS) \ + 0.5*(factor + logdetD) \ + 0.5*np.dot(np.transpose(self._j), mu) \ - - 0.5*np.dot(np.transpose(data), np.dot(np.diag(Wvalues), data)) + - 0.5*np.dot(np.transpose(np.conjugate(data)), np.dot(np.diag(Wvalues), data)) return -log_likelihood result = minimize(likelihood, x0=np.array([-5, 60, 1e5]), args=(V,), @@ -804,10 +808,9 @@ def inv_gamma_function(l, alpha, beta): method="Nelder-Mead", tol=1e-7, ) m, c, l = result.x - print("Result: ", "m: ", m, "c: ", c, "l: ", "{:e}".format(l)) + print("Best parameters: ", "m: ", m, "c: ", c, "l: ", "{:e}".format(l)) return [m, c, l] - def _fit(self): """Compute the mean and variance""" # Compute the inverse prior covariance, S(p)^-1 @@ -980,7 +983,8 @@ def num_fields(self): @property def size(self): """Number of points in reconstruction""" - return self._DHT.size + #return self._DHT.size + return self._DFT.size class LogNormalMAPModel: