diff --git a/CHANGELOG.md b/CHANGELOG.md index 8c5144ffbf..b35002e415 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,9 @@ All notable changes to this project will be documented in this file. ## [Unreleased] +### Added +- Support for MWA Average Embedded Element beams in UVBeam. + ## [3.2.6] - 2025-03-13 ### Added diff --git a/docs/uvbeam_tutorial.rst b/docs/uvbeam_tutorial.rst index 361924e586..035b95388e 100644 --- a/docs/uvbeam_tutorial.rst +++ b/docs/uvbeam_tutorial.rst @@ -69,9 +69,7 @@ UVBeam: Instantiating a UVBeam object from a file (i.e. reading data) Use the :meth:`pyuvdata.UVBeam.from_file` to instantiate a UVBeam object from data in a file (alternatively you can create an object with no inputs and then -call the :meth:`pyuvdata.UVBeam.read` method). Most file types require a single -file or folder to instantiate an object, FHD and raw MWA correlator data sets -require the user to specify multiple files for each dataset. +call the :meth:`pyuvdata.UVBeam.read` method). ``pyuvdata`` can also be used to create a UVBeam object from arrays in memory (see :ref:`new_uvbeam`) and to read in multiple datasets (files) into a single @@ -82,18 +80,29 @@ a) Instantiate an object from a single file or folder ***************************************************** BeamFITS and MWA beams are stored in a single file. -To get all the frequencies available for the MWA full embedded element beam -you need to download the output simulation file via -`wget http://cerberus.mwa128t.org/mwa_full_embedded_element_pattern.h5` -For this tutorial we use the file saved in the test data which only -contains a few frequencies. +To get all the frequencies available for the MWA beams you need to download the +output simulation file(s) for the model you want: + +- for the full embedded element (FEE) beam use: + ``wget http://cerberus.mwa128t.org/mwa_full_embedded_element_pattern.h5`` +- for the average embedded element (AEE) beam use: + ``wget https://github.com/MWATelescope/mwa_pb/blob/master/mwa_pb/data/Jmatrix.fits`` + and + ``wget https://github.com/MWATelescope/mwa_pb/blob/master/mwa_pb/data/Zmatrix.fits`` + +For this tutorial we use the files saved in the test data which only +contain a few frequencies. For MWA beams, which are composed of phased arrays of dipoles, you can specify delay and amplitude arrays to generate beams pointed any where or with varying gains per dipole. Set a delay to 32 to get a beam where that dipole is turned -off (or set the amplitude to zero). The native format of the beam is spherical -harmonic modes, so there is also an option `pixels_per_deg` to set the output -beam resolution (default is 5 pixels per degree). +off (or set the amplitude to zero). + +Depending on the beam version you want, you also need to specify some other paramers: + +- For the AEE beam you need to pass the Z matrix file to the ``mwa_zfile`` parameter. +- The native format of the FEE beam is spherical harmonic modes, so there is also + a ``pixels_per_deg`` parameter to set the output beam resolution (default is 5 pixels per degree). .. clear-namespace @@ -103,16 +112,20 @@ beam resolution (default is 5 pixels per degree). from pyuvdata import UVBeam from pyuvdata.datasets import fetch_data - filename = fetch_data("mwa_full_EE") + fee_filename = fetch_data("mwa_full_EE") + aee_filename = fetch_data("mwa_jmatrix") + aee_z_filename = fetch_data("mwa_zmatrix") # for a zenith pointed beam let the delays default to all zeros - beam = UVBeam.from_file(filename) + fee_beam = UVBeam.from_file(fee_filename) + aee_beam = UVBeam.from_file(aee_filename, mwa_zfile=aee_z_filename) # to point the beam off zenith apply a delay slope across the tile # use the same pointing for both pols. delays = np.empty((2, 16), dtype=int) for pol in range(2): delays[pol] = np.tile(np.arange(0, 8, 2), 4) - beam = UVBeam.from_file(filename, pixels_per_deg=1, delays=delays) + fee_beam = UVBeam.from_file(fee_filename, pixels_per_deg=1, delays=delays) + aee_beam = UVBeam.from_file(aee_filename, mwa_zfile=aee_z_filename, delays=delays) a) Instantiate an object from CST or FEKO beam files diff --git a/src/pyuvdata/data/mwa_config_data/mwa_lna_impedance.txt b/src/pyuvdata/data/mwa_config_data/mwa_lna_impedance.txt new file mode 100644 index 0000000000..0954a07253 --- /dev/null +++ b/src/pyuvdata/data/mwa_config_data/mwa_lna_impedance.txt @@ -0,0 +1,452 @@ +frequency_Hz Impedance_real Impedance_imag +50000000.0 52.803 151.76 +51000000.0 52.446 154.04 +52000000.0 52.19 157.03 +53000000.0 51.811 160.69 +54000000.0 51.528 164.85 +55000000.0 51.802 169.38 +56000000.0 52.157 174.17 +57000000.0 52.522 179.32 +58000000.0 53.377 184.66 +59000000.0 54.249 189.98 +60000000.0 55.512 195.65 +61000000.0 56.949 201.94 +62000000.0 58.519 208.08 +63000000.0 60.179 214.4 +64000000.0 62.755 220.82 +65000000.0 65.071 227.42 +66000000.0 67.949 234.58 +67000000.0 71.302 241.31 +68000000.0 74.795 248.26 +69000000.0 78.823 255.1 +70000000.0 83.465 262.45 +71000000.0 88.307 268.87 +72000000.0 93.864 275.2 +73000000.0 99.639 281.3 +74000000.0 105.58 287.5 +75000000.0 111.82 293.16 +76000000.0 118.57 298.85 +77000000.0 125.37 303.59 +78000000.0 132.08 308.64 +79000000.0 138.4 313.49 +80000000.0 144.44 318.07 +81000000.0 150.64 323.03 +82000000.0 157.06 327.82 +83000000.0 163.72 332.52 +84000000.0 170.03 337.71 +85000000.0 176.67 343.27 +86000000.0 184.35 349.36 +87000000.0 192.29 355.09 +88000000.0 200.02 360.97 +89000000.0 208.58 366.41 +90000000.0 216.59 372.96 +91000000.0 225.67 379.45 +92000000.0 235.53 386.17 +93000000.0 245.42 391.84 +94000000.0 256.16 398.14 +95000000.0 268.37 404.17 +96000000.0 281.23 410.06 +97000000.0 295.19 416.49 +98000000.0 308.93 420.98 +99000000.0 323.64 426.09 +100000000.0 339.04 431.74 +101000000.0 354.61 436.51 +102000000.0 371.38 440.57 +103000000.0 388.71 445.18 +104000000.0 406.16 447.12 +105000000.0 426.15 450.65 +106000000.0 446.62 452.77 +107000000.0 467.8 453.18 +108000000.0 489.21 450.95 +109000000.0 510.8 449.09 +110000000.0 534.27 446.18 +111000000.0 557.91 441.95 +112000000.0 582.01 436.09 +113000000.0 605.65 427.46 +114000000.0 629.14 418.13 +115000000.0 656.66 407.54 +116000000.0 681.79 394.13 +117000000.0 705.83 378.75 +118000000.0 730.36 361.53 +119000000.0 752.83 341.19 +120000000.0 775.31 320.56 +121000000.0 798.29 296.29 +122000000.0 816.81 272.67 +123000000.0 834.34 243.57 +124000000.0 851.51 217.39 +125000000.0 865.02 186.63 +126000000.0 876.49 157.17 +127000000.0 884.52 125.94 +128000000.0 891.0 96.191 +129000000.0 894.14 61.696 +130000000.0 897.08 27.864 +131000000.0 894.52 -3.5946 +132000000.0 892.67 -35.799 +133000000.0 887.2 -67.153 +134000000.0 881.6 -99.435 +135000000.0 872.62 -129.87 +136000000.0 864.69 -155.42 +137000000.0 853.69 -179.82 +138000000.0 840.15 -205.61 +139000000.0 825.47 -227.86 +140000000.0 811.32 -249.27 +141000000.0 794.84 -268.09 +142000000.0 778.56 -283.31 +143000000.0 760.98 -301.41 +144000000.0 745.14 -317.45 +145000000.0 728.93 -332.33 +146000000.0 712.62 -344.89 +147000000.0 695.47 -357.0 +148000000.0 678.85 -367.29 +149000000.0 662.72 -377.19 +150000000.0 645.2 -384.83 +151000000.0 628.67 -392.63 +152000000.0 611.52 -398.84 +153000000.0 594.87 -405.96 +154000000.0 577.85 -411.64 +155000000.0 561.2 -416.75 +156000000.0 545.01 -421.23 +157000000.0 531.35 -426.39 +158000000.0 516.58 -428.47 +159000000.0 501.43 -431.9 +160000000.0 486.88 -433.12 +161000000.0 473.71 -434.81 +162000000.0 461.38 -435.33 +163000000.0 448.13 -435.89 +164000000.0 434.9 -434.92 +165000000.0 421.63 -434.71 +166000000.0 409.75 -433.7 +167000000.0 397.58 -432.03 +168000000.0 386.65 -430.36 +169000000.0 376.25 -428.97 +170000000.0 366.39 -427.51 +171000000.0 356.49 -426.02 +172000000.0 347.71 -424.14 +173000000.0 339.13 -421.65 +174000000.0 330.76 -420.74 +175000000.0 322.54 -417.59 +176000000.0 314.44 -415.14 +177000000.0 305.26 -411.58 +178000000.0 297.84 -408.67 +179000000.0 290.45 -405.71 +180000000.0 282.98 -402.77 +181000000.0 276.15 -399.56 +182000000.0 268.95 -397.06 +183000000.0 262.05 -393.77 +184000000.0 256.48 -391.53 +185000000.0 250.55 -388.55 +186000000.0 244.67 -385.83 +187000000.0 239.24 -383.2 +188000000.0 234.1 -379.97 +189000000.0 229.17 -377.32 +190000000.0 224.07 -374.82 +191000000.0 219.45 -372.16 +192000000.0 214.22 -369.4 +193000000.0 209.86 -366.82 +194000000.0 205.51 -363.94 +195000000.0 201.02 -361.7 +196000000.0 196.64 -358.88 +197000000.0 192.06 -356.2 +198000000.0 187.99 -353.37 +199000000.0 184.01 -350.71 +200000000.0 180.82 -347.86 +201000000.0 177.31 -345.16 +202000000.0 173.86 -342.35 +203000000.0 170.04 -339.97 +204000000.0 167.37 -337.48 +205000000.0 163.77 -334.91 +206000000.0 160.71 -332.13 +207000000.0 157.16 -329.43 +208000000.0 153.89 -326.64 +209000000.0 150.9 -324.17 +210000000.0 148.42 -321.22 +211000000.0 145.32 -318.78 +212000000.0 142.8 -316.5 +213000000.0 140.51 -313.88 +214000000.0 138.07 -311.21 +215000000.0 135.65 -308.95 +216000000.0 133.56 -306.19 +217000000.0 131.42 -303.79 +218000000.0 129.42 -301.32 +219000000.0 127.57 -298.73 +220000000.0 125.73 -296.62 +221000000.0 123.7 -294.57 +222000000.0 121.82 -292.36 +223000000.0 119.93 -290.17 +224000000.0 118.27 -288.21 +225000000.0 116.39 -286.06 +226000000.0 114.48 -283.88 +227000000.0 112.76 -281.88 +228000000.0 111.13 -280.13 +229000000.0 109.49 -278.22 +230000000.0 107.74 -276.41 +231000000.0 105.99 -274.7 +232000000.0 104.74 -272.75 +233000000.0 103.35 -270.97 +234000000.0 101.76 -269.11 +235000000.0 100.21 -267.31 +236000000.0 98.957 -265.44 +237000000.0 97.601 -263.75 +238000000.0 96.383 -261.71 +239000000.0 94.926 -260.1 +240000000.0 93.655 -258.33 +241000000.0 92.223 -256.45 +242000000.0 91.028 -254.56 +243000000.0 89.743 -252.72 +244000000.0 88.719 -251.07 +245000000.0 87.393 -249.27 +246000000.0 86.092 -247.57 +247000000.0 84.946 -245.82 +248000000.0 83.73 -244.29 +249000000.0 82.535 -242.85 +250000000.0 81.653 -241.0 +251000000.0 80.563 -239.12 +252000000.0 79.475 -237.62 +253000000.0 78.791 -235.9 +254000000.0 77.636 -234.28 +255000000.0 76.866 -232.99 +256000000.0 76.39 -231.17 +257000000.0 75.255 -230.01 +258000000.0 74.317 -228.71 +259000000.0 73.493 -227.35 +260000000.0 72.697 -225.99 +261000000.0 72.096 -224.74 +262000000.0 71.661 -223.36 +263000000.0 70.576 -222.31 +264000000.0 69.961 -221.11 +265000000.0 69.4 -220.12 +266000000.0 68.796 -218.83 +267000000.0 67.94 -217.5 +268000000.0 67.31 -216.31 +269000000.0 66.216 -215.06 +270000000.0 65.584 -213.54 +271000000.0 64.877 -212.28 +272000000.0 63.893 -210.78 +273000000.0 63.306 -209.5 +274000000.0 62.675 -208.11 +275000000.0 61.585 -207.05 +276000000.0 61.274 -205.98 +277000000.0 60.893 -204.93 +278000000.0 60.162 -203.75 +279000000.0 59.491 -202.39 +280000000.0 58.785 -201.23 +281000000.0 58.037 -200.33 +282000000.0 57.763 -199.1 +283000000.0 57.082 -197.79 +284000000.0 56.345 -196.7 +285000000.0 55.75 -195.54 +286000000.0 55.49 -194.62 +287000000.0 54.995 -193.59 +288000000.0 54.358 -192.45 +289000000.0 53.946 -191.27 +290000000.0 53.429 -190.25 +291000000.0 52.941 -189.14 +292000000.0 52.774 -188.19 +293000000.0 52.502 -187.32 +294000000.0 52.019 -186.36 +295000000.0 51.827 -185.47 +296000000.0 51.392 -184.55 +297000000.0 50.93 -183.64 +298000000.0 50.507 -182.67 +299000000.0 50.048 -181.69 +300000000.0 49.391 -180.64 +301000000.0 48.997 -179.7 +302000000.0 48.551 -178.69 +303000000.0 48.098 -177.85 +304000000.0 47.737 -176.85 +305000000.0 47.302 -175.93 +306000000.0 46.857 -174.99 +307000000.0 46.454 -174.01 +308000000.0 46.187 -173.1 +309000000.0 45.907 -172.37 +310000000.0 45.59 -171.47 +311000000.0 45.215 -170.71 +312000000.0 44.848 -170.01 +313000000.0 44.296 -169.24 +314000000.0 43.968 -168.56 +315000000.0 43.492 -167.7 +316000000.0 43.203 -166.91 +317000000.0 42.812 -166.24 +318000000.0 42.453 -165.32 +319000000.0 42.256 -164.43 +320000000.0 42.245 -163.79 +321000000.0 41.847 -163.03 +322000000.0 41.537 -162.33 +323000000.0 41.211 -161.44 +324000000.0 40.928 -160.64 +325000000.0 40.545 -159.95 +326000000.0 40.022 -159.18 +327000000.0 39.594 -158.33 +328000000.0 39.436 -157.57 +329000000.0 39.083 -156.85 +330000000.0 38.585 -155.97 +331000000.0 38.308 -155.22 +332000000.0 38.104 -154.41 +333000000.0 37.944 -153.69 +334000000.0 37.686 -152.85 +335000000.0 37.312 -152.05 +336000000.0 37.102 -151.4 +337000000.0 37.074 -150.78 +338000000.0 36.776 -149.94 +339000000.0 36.544 -149.26 +340000000.0 36.359 -148.6 +341000000.0 36.08 -147.97 +342000000.0 35.99 -147.3 +343000000.0 35.746 -146.63 +344000000.0 35.507 -145.93 +345000000.0 35.38 -145.39 +346000000.0 35.17 -144.62 +347000000.0 34.934 -143.98 +348000000.0 34.799 -143.29 +349000000.0 34.627 -142.76 +350000000.0 34.442 -142.15 +351000000.0 34.097 -141.54 +352000000.0 33.835 -140.87 +353000000.0 33.512 -140.46 +354000000.0 33.297 -139.72 +355000000.0 32.991 -139.29 +356000000.0 32.688 -138.53 +357000000.0 32.479 -137.92 +358000000.0 32.298 -137.37 +359000000.0 32.054 -136.76 +360000000.0 31.97 -136.07 +361000000.0 31.669 -135.6 +362000000.0 31.519 -134.87 +363000000.0 31.281 -134.45 +364000000.0 31.115 -133.84 +365000000.0 30.888 -133.34 +366000000.0 30.746 -132.8 +367000000.0 30.581 -132.2 +368000000.0 30.414 -131.7 +369000000.0 30.129 -131.1 +370000000.0 30.007 -130.48 +371000000.0 29.826 -129.88 +372000000.0 29.702 -129.29 +373000000.0 29.507 -128.72 +374000000.0 29.274 -128.22 +375000000.0 29.115 -127.6 +376000000.0 28.968 -127.13 +377000000.0 28.758 -126.62 +378000000.0 28.543 -126.12 +379000000.0 28.386 -125.6 +380000000.0 28.211 -125.14 +381000000.0 28.065 -124.65 +382000000.0 27.857 -124.23 +383000000.0 27.758 -123.82 +384000000.0 27.619 -123.41 +385000000.0 27.564 -122.98 +386000000.0 27.532 -122.53 +387000000.0 27.377 -122.06 +388000000.0 27.215 -121.69 +389000000.0 27.072 -121.15 +390000000.0 26.882 -120.59 +391000000.0 26.655 -120.07 +392000000.0 26.442 -119.52 +393000000.0 26.172 -119.02 +394000000.0 25.979 -118.44 +395000000.0 25.805 -117.82 +396000000.0 25.732 -117.32 +397000000.0 25.542 -116.85 +398000000.0 25.435 -116.32 +399000000.0 25.251 -115.89 +400000000.0 24.992 -115.5 +401000000.0 24.86 -115.1 +402000000.0 24.668 -114.68 +403000000.0 24.544 -114.24 +404000000.0 24.343 -113.81 +405000000.0 24.136 -113.33 +406000000.0 23.995 -112.84 +407000000.0 23.9 -112.27 +408000000.0 23.732 -111.82 +409000000.0 23.654 -111.39 +410000000.0 23.413 -110.92 +411000000.0 23.39 -110.41 +412000000.0 23.214 -109.92 +413000000.0 23.097 -109.4 +414000000.0 22.958 -108.96 +415000000.0 22.836 -108.51 +416000000.0 22.7 -108.03 +417000000.0 22.618 -107.6 +418000000.0 22.428 -107.26 +419000000.0 22.419 -106.93 +420000000.0 22.344 -106.53 +421000000.0 22.267 -106.2 +422000000.0 22.137 -105.8 +423000000.0 22.05 -105.3 +424000000.0 21.98 -104.92 +425000000.0 21.851 -104.54 +426000000.0 21.762 -104.11 +427000000.0 21.648 -103.81 +428000000.0 21.493 -103.38 +429000000.0 21.433 -102.92 +430000000.0 21.34 -102.64 +431000000.0 21.154 -102.32 +432000000.0 20.983 -101.8 +433000000.0 20.861 -101.37 +434000000.0 20.606 -101.01 +435000000.0 20.523 -100.61 +436000000.0 20.314 -100.25 +437000000.0 20.112 -99.741 +438000000.0 19.995 -99.266 +439000000.0 19.955 -98.997 +440000000.0 19.758 -98.615 +441000000.0 19.688 -98.133 +442000000.0 19.582 -97.774 +443000000.0 19.581 -97.42 +444000000.0 19.475 -97.13 +445000000.0 19.387 -96.75 +446000000.0 19.285 -96.317 +447000000.0 19.257 -95.976 +448000000.0 19.19 -95.605 +449000000.0 19.042 -95.164 +450000000.0 18.918 -94.76 +451000000.0 18.813 -94.341 +452000000.0 18.675 -93.897 +453000000.0 18.573 -93.465 +454000000.0 18.45 -93.078 +455000000.0 18.416 -92.723 +456000000.0 18.373 -92.336 +457000000.0 18.222 -91.964 +458000000.0 18.151 -91.644 +459000000.0 18.142 -91.32 +460000000.0 18.12 -90.964 +461000000.0 18.066 -90.571 +462000000.0 17.921 -90.204 +463000000.0 17.831 -89.837 +464000000.0 17.799 -89.529 +465000000.0 17.705 -89.176 +466000000.0 17.555 -88.883 +467000000.0 17.417 -88.52 +468000000.0 17.354 -88.098 +469000000.0 17.344 -87.786 +470000000.0 17.25 -87.435 +471000000.0 17.144 -87.005 +472000000.0 17.013 -86.724 +473000000.0 16.952 -86.38 +474000000.0 16.872 -86.119 +475000000.0 16.786 -85.864 +476000000.0 16.636 -85.441 +477000000.0 16.569 -85.112 +478000000.0 16.519 -84.808 +479000000.0 16.461 -84.443 +480000000.0 16.327 -84.082 +481000000.0 16.251 -83.738 +482000000.0 16.066 -83.371 +483000000.0 15.943 -83.089 +484000000.0 15.814 -82.782 +485000000.0 15.669 -82.459 +486000000.0 15.522 -82.064 +487000000.0 15.458 -81.733 +488000000.0 15.293 -81.351 +489000000.0 15.243 -81.042 +490000000.0 15.146 -80.671 +491000000.0 15.05 -80.296 +492000000.0 14.973 -79.911 +493000000.0 14.955 -79.614 +494000000.0 14.876 -79.266 +495000000.0 14.839 -78.892 +496000000.0 14.723 -78.569 +497000000.0 14.683 -78.282 +498000000.0 14.643 -78.058 +499000000.0 14.614 -77.894 +500000000.0 14.594 -77.713 diff --git a/src/pyuvdata/data/test_data.yaml b/src/pyuvdata/data/test_data.yaml index d3724312a9..b123763841 100644 --- a/src/pyuvdata/data/test_data.yaml +++ b/src/pyuvdata/data/test_data.yaml @@ -43,6 +43,8 @@ mwa_cotter_phase_test2: visibility_data/MWA/1133866760_rephase.uvfits mwa_fhd: visibility_data/MWA/fhd_vis_data.tar.gz mwa_fhd_cal: calibration_solutions/MWA/fhd_cal_data.tar.gz mwa_full_EE: primary_beams/MWA/mwa_full_EE_test.h5 +mwa_jmatrix: primary_beams/MWA/JMatrix_3freq.fits +mwa_zmatrix: primary_beams/MWA/ZMatrix_3freq.fits mwax_2021_metafits: visibility_data/MWA/mwa_corr_fits_testfiles/1320409688.metafits mwax_2021_raw_gpubox: visibility_data/MWA/mwa_corr_fits_testfiles/1320409688_20211108122750_mini_ch137_000.fits ovro_lwa: visibility_data/LWA/2018-03-21-01_26_33_0004384620257280_000000_downselected.ms.tar.gz diff --git a/src/pyuvdata/data/test_data_registry.txt b/src/pyuvdata/data/test_data_registry.txt index eda7b29cf0..d7b73c8e0f 100644 --- a/src/pyuvdata/data/test_data_registry.txt +++ b/src/pyuvdata/data/test_data_registry.txt @@ -68,4 +68,6 @@ primary_beams/HERA/HERABEAM.FITS a28232d8f63af528afae7ca529856809d478db86b3b3b7d primary_beams/LWA/OVRO_LWA_x.ffe c69e5ac1afa638d2a80663ded0bf4d69c1273a55bc16b52e15e3d96209a3087e primary_beams/LWA/OVRO_LWA_y.ffe 3e2f165ac82c7f422815c7b9911ae61dc7ed43dff00556ffd59b58248a93119a primary_beams/MWA/mwa_full_EE_test.h5 9246b9e98dc99d4284e7865ebae221b43aaef64893830cad00587d265eef2bfe +primary_beams/MWA/JMatrix_3freq.fits c1a0100a3c0a6999d32c393f0712cb9d86f6f4340376d4cfea9eba120c94b4c0 +primary_beams/MWA/ZMatrix_3freq.fits fe79c5601cab0d4cd82260860ad98f66e0dcf7098b51c5b51108d344a76654a6 flagsets/HERA/zen.2457698.40355.xx.HH.uvcAA.testuvflag.h5 20ff14975bcedf5edd447e3c5db3456d9298846d7e5da5bb084c7b0662327875 diff --git a/src/pyuvdata/datasets.py b/src/pyuvdata/datasets.py index 5208290941..faf789c2ec 100644 --- a/src/pyuvdata/datasets.py +++ b/src/pyuvdata/datasets.py @@ -12,7 +12,7 @@ from pyuvdata.data import DATA_PATH # set the data version. This should be updated when the test data are changed. -data_version = "v0.0.4" +data_version = "v0.0.5" pup = pooch.create( # Use the default cache folder for the operating system diff --git a/src/pyuvdata/uvbeam/mwa_beam.py b/src/pyuvdata/uvbeam/mwa_beam.py index 10112f14e7..ecebbd4dbd 100644 --- a/src/pyuvdata/uvbeam/mwa_beam.py +++ b/src/pyuvdata/uvbeam/mwa_beam.py @@ -7,15 +7,27 @@ import h5py import numpy as np +from astropy import constants +from astropy.io import fits from docstring_parser import DocstringStyle from scipy.special import factorial, lpmv # associated Legendre function from .. import utils +from ..data import DATA_PATH from ..docstrings import copy_replace_short_description from . import UVBeam __all__ = ["P1sin", "P1sin_array", "MWABeam"] +# dipole spacing in meters +MWA_DIPOLE_SPACING_M = 1.1 + +# 435 picoseconds is base delay length unit +MWA_BASE_DELAY_S = 4.35e-10 + +MWA_NFEED = 2 +MWA_NDIPOLE = 16 + def P1sin(nmax, theta): # noqa N802 """ @@ -192,6 +204,44 @@ def P1sin_array(nmax, theta): # noqa N802 return P_sin.transpose(), P1.transpose() +def _get_freq_inds_use(freqs_hz, freq_range=None): + """ + Get the frequencies to use. + + Parameters + ---------- + freqs_hz : array of float + Frequencies in file. + freq_range : tuple of float in Hz + If given, the lower and upper limit of the frequencies to read in. Default + is to use all frequencies. + + Returns + ------- + freq_inds_use : array of int + Indices into the freqs_hz to use. + """ + if freq_range is not None: + if np.array(freq_range).size != 2: + raise ValueError("freq_range must have 2 elements.") + freqs_mask = utils.tools._is_between( + freqs_hz, np.asarray(freq_range)[np.newaxis] + ) + freq_inds_use = np.nonzero(freqs_mask)[0] + if freq_inds_use.size < 1: + raise ValueError( + "No frequencies available in freq_range. " + f"Available frequencies are between {np.min(freqs_hz)} Hz " + f"and {np.max(freqs_hz)} Hz" + ) + if freq_inds_use.size < 2: + warnings.warn("Only one available frequency in freq_range.") + else: + freq_inds_use = np.arange(freqs_hz.size) + + return freq_inds_use + + class MWABeam(UVBeam): """ Defines an MWA-specific subclass of UVBeam for representing MWA beams. @@ -199,8 +249,11 @@ class MWABeam(UVBeam): This class should not be interacted with directly, instead use the read_mwa_beam method on the UVBeam class. - This is based on https://github.com/MWATelescope/mwa_pb/ but we don't import - that module because it's not python 3 compatible. + This is based on https://github.com/MWATelescope/mwa_pb/ and offers support + for either the Fully Embedded Element (FEE) model or the Average Embedded + Element (AEE) model. The FEE model is the most current model and is generally + considered to be the best match to the instrument beam, but we provide the + older AEE model for comparison as well. Note that the azimuth convention for the UVBeam object is different than the azimuth convention in the mwa_pb repo. In that repo, the azimuth convention @@ -227,7 +280,7 @@ def _read_metadata(self, h5filepath): ------- freqs_hz : array of int Frequencies in Hz present in the file. - pol_names : list of str + feed_names : list of str Polarizations present in the file. dipole_names : Dipoles names present in the file. @@ -235,7 +288,7 @@ def _read_metadata(self, h5filepath): Dictionary keyed on pol and freq, giving max number of modes in the file for each pol and freq. """ - pol_names = set() + feed_names = set() dipole_names = set() freqs_hz = set() other_names = [] @@ -245,7 +298,7 @@ def _read_metadata(self, h5filepath): if name.startswith("X") or name.startswith("Y"): pol = name[0] dipole, freq = name[1:].split("_") - pol_names.add(pol) + feed_names.add(pol) dipole_names.add(dipole) freq = np.int64(freq) freqs_hz.add(freq) @@ -263,20 +316,20 @@ def _read_metadata(self, h5filepath): else: other_names.append(name) - pol_names = sorted(pol_names) + feed_names = sorted(feed_names) dipole_names = np.asarray(sorted(dipole_names, key=int)) freqs_hz = np.array(sorted(freqs_hz)) - return freqs_hz, pol_names, dipole_names, max_length + return freqs_hz, feed_names, dipole_names, max_length def _get_beam_modes( self, *, h5filepath, freqs_hz, - pol_names, + feed_names, dipole_names, max_length, delays, @@ -292,7 +345,7 @@ def _get_beam_modes( harmonic modes. freqs_hz : array of int Frequencies in Hz to get modes for. Must be present in the file. - pol_names : list of str + feed_names : list of str Polarizations to get modes for. Must be present in the file. dipole_names : array of str Dipoles names present in the file. @@ -300,10 +353,10 @@ def _get_beam_modes( Dictionary keyed on pol and freq, giving max number of modes in the file for each pol and freq. delays : array of ints - Array of MWA beamformer delay steps. Should be shape (n_pols, n_dipoles). + Array of MWA beamformer delay steps. Should be shape (n_feeds, n_dipoles). amplitudes : array of floats Array of dipole amplitudes, these are absolute values - (i.e. relatable to physical units). Should be shape (n_pols, n_dipoles). + (i.e. relatable to physical units). Should be shape (n_feeds, n_dipoles). Returns ------- @@ -311,7 +364,7 @@ def _get_beam_modes( A multi-level dict keyed on (in order) pol, freq, mode name (Q1, Q2, M, N). """ beam_modes = {} - for pol_i, pol in enumerate(pol_names): + for pol_i, pol in enumerate(feed_names): beam_modes[pol] = {} for freq in freqs_hz: # Calculate complex excitation voltages @@ -386,7 +439,7 @@ def _get_beam_modes( } return beam_modes - def _get_response(self, *, freqs_hz, pol_names, beam_modes, phi_arr, theta_arr): + def _get_response(self, *, freqs_hz, feed_names, beam_modes, az_array, za_array): """ Calculate full Jones matrix response (E-field) of beam on a regular az/za grid. @@ -394,13 +447,13 @@ def _get_response(self, *, freqs_hz, pol_names, beam_modes, phi_arr, theta_arr): ---------- freqs_hz : array of int Frequencies in Hz to get modes for. Must be present in the file. - pol_names : list of str + feed_names : list of str Polarizations to get modes for. Must be present in the file. beam_modes : dict A multi-level dict keyed on (in order) pol, freq, mode name (Q1, Q2, M, N). - phi_arr : float or array of float + az_array : float or array of float azimuth angles (radians), east through north. - theta_arr : float or array of float + za_array : float or array of float zenith angles (radian) Returns @@ -412,11 +465,11 @@ def _get_response(self, *, freqs_hz, pol_names, beam_modes, phi_arr, theta_arr): """ jones = np.zeros( - (len(pol_names), 2, freqs_hz.size, phi_arr.size, theta_arr.size), + (len(feed_names), 2, freqs_hz.size, az_array.size, za_array.size), dtype=np.complex128, ) - for pol_i, pol in enumerate(pol_names): + for pol_i, pol in enumerate(feed_names): for freq_i, freq in enumerate(freqs_hz): M = beam_modes[pol][freq]["M"] N = beam_modes[pol][freq]["N"] @@ -451,11 +504,11 @@ def _get_response(self, *, freqs_hz, pol_names, beam_modes, phi_arr, theta_arr): # theta and phi are direction coordinates phi_comp = np.ascontiguousarray( - np.exp(1.0j * np.outer(phi_arr, range(-nmax, nmax + 1))) + np.exp(1.0j * np.outer(az_array, range(-nmax, nmax + 1))) ) - (P_sin, P1) = P1sin_array(nmax, theta_arr) - M_u = np.outer(np.cos(theta_arr), np.abs(M)) + (P_sin, P1) = P1sin_array(nmax, za_array) + M_u = np.outer(np.cos(za_array), np.abs(M)) phi_const = C_MN * MabsM / (N * (N + 1)) ** 0.5 emn_T = ( @@ -469,12 +522,8 @@ def _get_response(self, *, freqs_hz, pol_names, beam_modes, phi_arr, theta_arr): # Use a matrix multiplication to calculate Emn_P and Emn_T. # Sum results of Emn_P and emn_T for each unique M - emn_P_sum = np.zeros( - (len(theta_arr), 2 * nmax + 1), dtype=np.complex128 - ) - emn_T_sum = np.zeros( - (len(theta_arr), 2 * nmax + 1), dtype=np.complex128 - ) + emn_P_sum = np.zeros((len(za_array), 2 * nmax + 1), dtype=np.complex128) + emn_T_sum = np.zeros((len(za_array), 2 * nmax + 1), dtype=np.complex128) for m in range(-nmax, nmax + 1): emn_P_sum[:, m + nmax] = np.sum(emn_P[:, m == M], axis=1) emn_T_sum[:, m + nmax] = np.sum(emn_T[:, m == M], axis=1) @@ -490,8 +539,7 @@ def _get_response(self, *, freqs_hz, pol_names, beam_modes, phi_arr, theta_arr): return jones - @copy_replace_short_description(UVBeam.read_mwa_beam, style=DocstringStyle.NUMPYDOC) - def read_mwa_beam( + def _read_fee_jones( self, h5filepath, *, @@ -499,11 +547,6 @@ def read_mwa_beam( amplitudes=None, pixels_per_deg=5, freq_range=None, - run_check=True, - check_extra=True, - run_check_acceptability=True, - check_auto_power=True, - fix_auto_power=True, ): """Read in the full embedded element MWA beam.""" # update filename attribute @@ -511,25 +554,325 @@ def read_mwa_beam( self.filename = [basename] self._filename.form = (1,) - freqs_hz, pol_names, dipole_names, max_length = self._read_metadata(h5filepath) + freqs_hz, feed_names, dipole_names, max_length = self._read_metadata(h5filepath) + + freqs_inds_use = _get_freq_inds_use(freqs_hz, freq_range=freq_range) + freqs_use = freqs_hz[freqs_inds_use] + + beam_modes = self._get_beam_modes( + h5filepath=h5filepath, + freqs_hz=freqs_hz, + feed_names=feed_names, + dipole_names=dipole_names, + max_length=max_length, + delays=delays, + amplitudes=amplitudes, + ) + + n_phi = np.floor(360 * pixels_per_deg) + n_theta = np.floor(90 * pixels_per_deg) + 1 + za_array = np.deg2rad(np.arange(0, n_theta) / pixels_per_deg) + az_array = np.deg2rad(np.arange(0, n_phi) / pixels_per_deg) + + jones = self._get_response( + freqs_hz=freqs_use, + feed_names=feed_names, + beam_modes=beam_modes, + az_array=az_array, + za_array=za_array, + ) + + # The array that come from `_get_response` has shape shape + # (Npol, 2, Nfreq, Nphi, Ntheta) + # UVBeam wants shape + # ('Naxes_vec', 1, 'Nfeeds', 'Nfreqs', 'Naxes2', 'Naxes1') + # where the Naxes_vec dimension lines up with the 2 from `_get_response`, + # Nfeeds is UVBeam's Npol for E-field beams, + # and axes (2, 1) correspond to (theta, phi) + # Then add an empty dimension for Nspws. + jones = np.transpose(jones, axes=[1, 0, 2, 4, 3]) + + return jones, freqs_use, za_array, az_array + + def _read_aee_jones( + self, + jonesfile, + *, + zfile, + delays=None, + amplitudes=None, + freq_range=None, + include_cross_feed_coupling=True, + ): + """Read in the average embedded element MWA beam.""" + # update filename attribute + jones_basename = os.path.basename(jonesfile) + z_basename = os.path.basename(zfile) + self.filename = [jones_basename, z_basename] + self._filename.form = (2,) + + n_feed = MWA_NFEED + n_dipole = MWA_NDIPOLE + + with fits.open(jonesfile) as jfile: + n_freqs = len(jfile) + freqs_hz = np.ndarray(n_freqs, dtype=float) + for f_ind in range(n_freqs): + freqs_hz[f_ind] = jfile[f_ind].header["freq"] + raw_theta = jfile[0].data[:, 0] + raw_phi = jfile[0].data[:, 1] + + freqs_inds_use = _get_freq_inds_use(freqs_hz, freq_range=freq_range) + freqs_use = freqs_hz[freqs_inds_use] + n_freqs = freqs_inds_use.size + + theta = np.unique(raw_theta) + phi = np.unique(raw_phi) + + n_theta = theta.size + n_phi = phi.size + if not n_theta * n_phi == raw_theta.size: + raise ValueError("Data does not appear to be on a grid") + + phi_grid, theta_grid = np.meshgrid(phi, theta) + + if not np.allclose(raw_theta.reshape(n_phi, n_theta).T, theta_grid): + raise ValueError("thetas do not appear to be on expected grid") + if not np.allclose(raw_phi.reshape(n_phi, n_theta).T, phi_grid): + raise ValueError("phis do not appear to be on expected grid") + + # convert theta, phi to radians, rename + az_grid = np.deg2rad(phi_grid) + za_grid = np.deg2rad(theta_grid) + az_array = np.deg2rad(phi) + za_array = np.deg2rad(theta) + + # this is just the dipole jones for now, but will be updated later + aee_jones = np.ndarray((2, 2, n_freqs, n_theta, n_phi), dtype=complex) + + with fits.open(jonesfile) as jfile: + # This is the comment from the JMatrix FITS file: + # Cols: theta phi real(Jxt(t,p)) imag(Jxt(t,p)) real(Jxp(t,p)) + # imag(Jxp(t,p)) real(Jyt(t,p)) imag(Jyt(t,p)) real(Jyp(t,p)) + # imag(Jyp(t,p))) + # Where theta is the zenith angle, phi is angle measured clockwise + # from +east direction looking down + # Jxt is the Jones mapping unit vec in theta (t) direction to the x + # (east-west) dipole etc + + for fi_arr, f_ind in enumerate(freqs_inds_use): + data = jfile[f_ind].data + + if not np.allclose(raw_theta, data[:, 0]): + raise ValueError("Inconsistent theta values across frequencies") + if not np.allclose(raw_phi, data[:, 1]): + raise ValueError("Inconsistent phi values across frequencies") + + aee_jones[1, 0, fi_arr] = ( + (data[:, 2] + 1j * data[:, 3]).reshape(n_phi, n_theta).T + ) + aee_jones[0, 0, fi_arr] = ( + (data[:, 4] + 1j * data[:, 5]).reshape(n_phi, n_theta).T + ) + aee_jones[1, 1, fi_arr] = ( + (data[:, 6] + 1j * data[:, 7]).reshape(n_phi, n_theta).T + ) + aee_jones[0, 1, fi_arr] = ( + (data[:, 8] + 1j * data[:, 9]).reshape(n_phi, n_theta).T + ) + + # Now we have to work out the apperature array factor -- the factor + # to multiply the average dipole jones by to account for all the + # apperature array stuff: geometric dipole delays, delay line settings, + # dipole impedances, dipole couplings + + # set up the dipole centers on a 4x4 grid with the dipole spacing for + # geometric delay calculation + x_centers, y_centers = np.meshgrid( + np.arange(4) * MWA_DIPOLE_SPACING_M, + np.flip(np.arange(4)) * MWA_DIPOLE_SPACING_M, + ) + x_centers = (x_centers - x_centers.mean()).flatten() + y_centers = (y_centers - y_centers.mean()).flatten() + z_centers = np.zeros(16, dtype=float) + + # calculate delays for dipoles relative to the center of the tile + # az runs from East to North + east_direction_cos = np.sin(za_grid) * np.cos(az_grid) + north_direction_cos = np.sin(za_grid) * np.sin(az_grid) + radial_direction_cos = np.cos(za_grid) + + geometric_dipole_delay = ( + np.outer(east_direction_cos.flatten(), x_centers) + + np.outer(north_direction_cos.flatten(), y_centers) + + np.outer(radial_direction_cos.flatten(), z_centers) + ) + + delays_sec = delays * MWA_BASE_DELAY_S + + lna_impedance_file = os.path.join( + DATA_PATH, "mwa_config_data", "mwa_lna_impedance.txt" + ) + + # the LNA impedances were measured on a different grid than the beam file + # interpolate the impedances to the frequencies we care about + data = np.genfromtxt(lna_impedance_file, skip_header=1) + impedance_freq = data[:, 0] + lna_impedance = data[:, 1] + 1j * data[:, 2] + lna_impedance_use = np.interp( + freqs_use, impedance_freq, np.real(lna_impedance) + ) + 1j * np.interp(freqs_use, impedance_freq, np.imag(lna_impedance)) + + # read the dipole coupling out of the Zmatrix file. these are 32 x 32 + # as all dipole pols can couple to all other dipole pols + dipole_coupling = np.zeros((n_freqs, n_dipole * 2, n_dipole * 2), dtype=complex) + with fits.open(zfile) as zf: + # Comment from ZMatrix FITS file and related code: + # Data are 32x32x2 cubes per ext. 1st plane Mag, 2nd plane phase + # ordering in Z matrix is 0-15:Y, 16-31:X + if not len(zf) == freqs_hz.size: + raise ValueError( + "Zmatrix file does not have as the same number of frequencies " + "as Jmatrix file." + ) + for fi_arr, f_ind in enumerate(freqs_inds_use): + if not np.isclose( + zf[f_ind].header["freq"], + freqs_hz[f_ind], + rtol=self._freq_array.tols[0], + atol=self._freq_array.tols[1], + ): + raise ValueError( + f"Zmatrix {f_ind}th freq does not match Jmatrix file." + ) + data = zf[f_ind].data + # convert to a proper complex number + data = data[0, :, :] * ( + np.cos(data[1, :, :]) + 1j * np.sin(data[1, :, :]) + ) + # move blocks around so we have X dipoles then Y dipoles + dipole_coupling[fi_arr, :n_dipole, :n_dipole] = data[ + n_dipole:, n_dipole: + ] + dipole_coupling[fi_arr, :n_dipole, n_dipole:] = data[ + n_dipole:, :n_dipole + ] + dipole_coupling[fi_arr, n_dipole:, :n_dipole] = data[ + :n_dipole, n_dipole: + ] + dipole_coupling[fi_arr, n_dipole:, n_dipole:] = data[ + :n_dipole, :n_dipole + ] + + n_couple = n_dipole * 2 + if not include_cross_feed_coupling: + # this is what FHD does. But by construction it excludes x <-> y coupling + # because rather than a 32 x 32 matrix we have (2 x 16 x 16) + n_couple = n_dipole + same_feed_coupling = np.zeros( + (2, n_freqs, n_dipole, n_dipole), dtype=complex + ) + same_feed_coupling[0] = dipole_coupling[:, :n_dipole, :n_dipole] + same_feed_coupling[1] = dipole_coupling[:, n_dipole:, n_dipole:] + dipole_coupling = same_feed_coupling + + z_total = dipole_coupling.copy() + for fi in range(n_freqs): + if include_cross_feed_coupling: + z_total[fi] += np.eye(n_couple) * lna_impedance_use[fi] + else: + for feed_i in range(n_feed): + z_total[feed_i, [fi]] += np.eye(n_couple) * lna_impedance_use[fi] + + # invert z to calculate currents (as voltage/impedance) + inv_z = np.linalg.inv(z_total) + + # signs in delayline term were adjusted to make the pointings go the right way + apparray_factor = np.zeros((2, n_freqs, n_theta, n_phi), dtype=complex) + for fi, this_f in enumerate(freqs_use): + k_conv = (2.0 * np.pi) * (this_f / constants.c.to("m/s").value) + delayline_term = np.exp( + -1j * 2.0 * np.pi * delays_sec * this_f * amplitudes + ) + if include_cross_feed_coupling: + port_current = np.dot(inv_z[fi], delayline_term.reshape(32)).reshape( + 2, 16 + ) + + for feed_i in range(n_feed): + if include_cross_feed_coupling: + port_current_use = port_current[feed_i] + else: + port_current_use = np.dot(inv_z[feed_i, fi], delayline_term[feed_i]) + + for dip_i in range(n_dipole): + geometric_dipole_factor = np.exp( + 1j * k_conv * geometric_dipole_delay[:, dip_i] + ) + # sum product of geometric factor & port current across dipoles + apparray_factor[feed_i, fi] += ( + geometric_dipole_factor * port_current_use[dip_i] + ).reshape((n_theta, n_phi)) + + for vec_i in range(2): + aee_jones[vec_i] *= apparray_factor + + return aee_jones, freqs_use, za_array, az_array + + @copy_replace_short_description(UVBeam.read_mwa_beam, style=DocstringStyle.NUMPYDOC) + def read_mwa_beam( + self, + filename, + *, + model_type=None, + zfile=None, + delays=None, + amplitudes=None, + pixels_per_deg=5, + freq_range=None, + include_cross_feed_coupling=True, + run_check=True, + check_extra=True, + run_check_acceptability=True, + check_auto_power=True, + fix_auto_power=True, + ): + """Read in MWA beam.""" + n_feed = MWA_NFEED + n_dipole = MWA_NDIPOLE + + if model_type is None: + _, extension = os.path.splitext(filename) + if extension == ".fits": + model_type = "aee" + elif extension == ".hdf5" or extension == ".h5": + model_type = "fee" + else: + raise ValueError( + "model_type could not be determined for MWA beam file, use " + "the model_type keyword to specify the type (one of 'fee' " + f"or 'aee'). Filename is: {filename}" + ) - n_dp = dipole_names.size - n_pol = len(pol_names) + if model_type == "aee" and zfile is None: + raise ValueError( + "zfile must be supplied for average embedded element MWA beam." + ) if delays is None: - delays = np.zeros([n_pol, n_dp], dtype="int") + delays = np.zeros([n_feed, n_dipole], dtype="int") else: if not np.issubdtype(delays.dtype, np.integer): raise ValueError("Delays must be integers.") if amplitudes is None: - amplitudes = np.ones([n_pol, n_dp]) + amplitudes = np.ones([n_feed, n_dipole]) - if amplitudes.shape != (n_pol, n_dp): - raise ValueError(f"amplitudes must be shape ({n_pol}, {n_dp})") + if amplitudes.shape != (n_feed, n_dipole): + raise ValueError(f"amplitudes must be shape ({n_feed}, {n_dipole})") - if delays.shape != (n_pol, n_dp): - raise ValueError(f"delays must be shape ({n_pol}, {n_dp})") + if delays.shape != (n_feed, n_dipole): + raise ValueError(f"delays must be shape ({n_feed}, {n_dipole})") if (delays > 32).any(): raise ValueError(f"There are delays greater than 32: {delays}") @@ -545,54 +888,31 @@ def read_mwa_beam( delays[terminated] = 0 amplitudes[terminated] = 0 - if freq_range is not None: - if np.array(freq_range).size != 2: - raise ValueError("freq_range must have 2 elements.") - freqs_use = freqs_hz[ - np.nonzero((freqs_hz >= freq_range[0]) & (freqs_hz <= freq_range[1])) - ] - if freqs_use.size < 1: - raise ValueError( - "No frequencies available in freq_range. " - f"Available frequencies are between {np.min(freqs_hz)} Hz " - f"and {np.max(freqs_hz)} Hz" - ) - if freqs_use.size < 2: - warnings.warn("Only one available frequency in freq_range.") + if model_type == "fee": + jones, freq_array, za_array, az_array = self._read_fee_jones( + filename, + delays=delays, + amplitudes=amplitudes, + pixels_per_deg=pixels_per_deg, + freq_range=freq_range, + ) else: - freqs_use = freqs_hz - - beam_modes = self._get_beam_modes( - h5filepath=h5filepath, - freqs_hz=freqs_hz, - pol_names=pol_names, - dipole_names=dipole_names, - max_length=max_length, - delays=delays, - amplitudes=amplitudes, - ) - - n_phi = np.floor(360 * pixels_per_deg) - n_theta = np.floor(90 * pixels_per_deg) + 1 - theta_arr = np.deg2rad(np.arange(0, n_theta) / pixels_per_deg) - phi_arr = np.deg2rad(np.arange(0, n_phi) / pixels_per_deg) - - jones = self._get_response( - freqs_hz=freqs_use, - pol_names=pol_names, - beam_modes=beam_modes, - phi_arr=phi_arr, - theta_arr=theta_arr, - ) - - # work out zenith normalization - # (MWA beams are peak normalized to 1 when pointed at zenith) - + jones, freq_array, za_array, az_array = self._read_aee_jones( + filename, + zfile=zfile, + delays=delays, + amplitudes=amplitudes, + freq_range=freq_range, + include_cross_feed_coupling=include_cross_feed_coupling, + ) # start filling in the object self.telescope_name = "MWA" self.feed_name = "MWA" self.feed_version = "1.0" - self.model_name = "full embedded element" + if model_type == "fee": + self.model_name = "full embedded element" + else: + self.model_name = "average embedded element" self.model_version = "1.0" self.history = ( "Sujito et al. full embedded element beam, derived from " @@ -601,7 +921,7 @@ def read_mwa_beam( delay_str_list = [] gain_str_list = [] - for pol in range(n_pol): + for pol in range(n_feed): delay_str_list.append( "[" + ", ".join([str(x) for x in delays[pol, :]]) + "]" ) @@ -620,7 +940,7 @@ def read_mwa_beam( self._set_efield() self.Naxes_vec = 2 self.Ncomponents_vec = 2 - self.feed_array = np.array([str(pol.lower()) for pol in pol_names]) + self.feed_array = np.array(["x", "y"]) self.feed_angle = np.where(self.feed_array == "x", np.pi / 2, 0.0) self.Nfeeds = self.feed_array.size self.mount_type = "phased" @@ -631,27 +951,19 @@ def read_mwa_beam( # to make the beam self.antenna_type = "simple" - self.Nfreqs = freqs_use.size - self.freq_array = freqs_use.astype(np.float64) + self.Nfreqs = freq_array.size + self.freq_array = freq_array.astype(np.float64) self.bandpass_array = np.ones(self.Nfreqs) self.pixel_coordinate_system = "az_za" self._set_cs_params() - self.axis1_array = phi_arr + self.axis1_array = az_array self.Naxes1 = self.axis1_array.size - self.axis2_array = theta_arr + self.axis2_array = za_array self.Naxes2 = self.axis2_array.size - # The array that come from `_get_response` has shape shape - # (Npol, 2, Nfreq, Nphi, Ntheta) - # UVBeam wants shape - # ('Naxes_vec', 1, 'Nfeeds', 'Nfreqs', 'Naxes2', 'Naxes1') - # where the Naxes_vec dimension lines up with the 2 from `_get_response`, - # Nfeeds is UVBeam's Npol for E-field beams, - # and axes (2, 1) correspond to (theta, phi) - # Then add an empty dimension for Nspws. - self.data_array = np.transpose(jones, axes=[1, 0, 2, 4, 3]) + self.data_array = jones self.basis_vector_array = np.zeros( (self.Naxes_vec, self.Ncomponents_vec, self.Naxes2, self.Naxes1) diff --git a/src/pyuvdata/uvbeam/uvbeam.py b/src/pyuvdata/uvbeam/uvbeam.py index 5728f90aba..3940e69580 100644 --- a/src/pyuvdata/uvbeam/uvbeam.py +++ b/src/pyuvdata/uvbeam/uvbeam.py @@ -4226,9 +4226,27 @@ def read_feko_beam(self, filename, **kwargs): self._convert_from_filetype(feko_beam_obj) del feko_beam_obj - def read_mwa_beam(self, h5filepath, **kwargs): + def read_mwa_beam(self, filename, **kwargs): """ - Read in the full embedded element MWA beam. + Read in either the full or averaged embedded element MWA beam. + + This is based on https://github.com/MWATelescope/mwa_pb/ and offers support + for either the Full Embedded Element (FEE) model or the Average Embedded + Element (AEE) model. The FEE model is the most current model and is generally + considered to be the best match to the instrument beam, but we provide the + older AEE model for comparison as well. + + For the FEE model, pass the h5 file containing the MWA full embedded + element spherical harmonic modes. Download via + `wget http://ws.mwatelescope.org/static/mwa_full_embedded_element_pattern.h5` + + For the AEE model, you need both the JMatrix.fits file and the ZMatrix.fits + files. Download via + `wget https://github.com/MWATelescope/mwa_pb/blob/master/mwa_pb/data/Jmatrix.fits` + and + `wget https://github.com/MWATelescope/mwa_pb/blob/master/mwa_pb/data/Zmatrix.fits` + Pass the JMatrix.fits file to the filename parameter and the ZMatrix.fits + file to the zfile parameter. Note that the azimuth convention in for the UVBeam object is different than the azimuth convention in the mwa_pb repo. In that repo, the azimuth convention is @@ -4238,11 +4256,16 @@ def read_mwa_beam(self, h5filepath, **kwargs): Parameters ---------- - h5filepath : str - path to input h5 file containing the MWA full embedded element spherical - harmonic modes. Download via - `wget http://ws.mwatelescope.org/static/mwa_full_embedded_element_pattern.h5` - (This reader is based on https://github.com/MWATelescope/mwa_pb). + filename : str + Path to the beam file containing the jones matrix information, either + the mwa_full_embedded_element_pattern.h5 for the FEE model or the + JMatrix.fits for the AEE model. + model_type : str + Model type to use, one of 'fee' or 'aee'. Will be inferred from + filename if not set and if possible. + zfile : str + Path to the file containing the antenna coupling matrix (ZMatrix.fits), + only used for the AEE model. delays : array of ints Array of MWA beamformer delay steps. Should be shape (n_pols, n_dipoles). amplitudes : array of floats @@ -4250,10 +4273,16 @@ def read_mwa_beam(self, h5filepath, **kwargs): (i.e. relatable to physical units). Should be shape (n_pols, n_dipoles). pixels_per_deg : float - Number of theta/phi pixels per degree. Sets the resolution of the beam. + Number of theta/phi pixels per degree. Sets the resolution of the beam, + only used for the FEE model. freq_range : array_like of float Range of frequencies to include in Hz, defaults to all available frequencies. Must be length 2. + include_cross_feed_coupling : bool + Option to include the couplings between the different feed directions + (i.e. couplings between x and y feeds). Including these couplings + is more correct, but the option is supplied to enable matching with + some codes that exclude these couplings. Default is True. run_check : bool Option to check for the existence and proper shapes of required parameters after reading in the file. @@ -4273,7 +4302,7 @@ def read_mwa_beam(self, h5filepath, **kwargs): from . import mwa_beam mwabeam_obj = mwa_beam.MWABeam() - mwabeam_obj.read_mwa_beam(h5filepath, **kwargs) + mwabeam_obj.read_mwa_beam(filename, **kwargs) self._convert_from_filetype(mwabeam_obj) del mwabeam_obj @@ -4314,9 +4343,12 @@ def read( extra_keywords=None, frequency_select=None, # mwa parameters + mwa_model_type=None, + mwa_zfile=None, delays=None, amplitudes=None, pixels_per_deg=5, + mwa_include_cross_feed_coupling=True, ): """ Read a generic file into a UVBeam object. @@ -4537,6 +4569,12 @@ def read( MWA --- + mwa_model_type : str + Model type to use, one of 'fee' or 'aee'. Will be inferred from + filename if not set and if possible. + mwa_zfile : str + Path to the file containing the antenna coupling matrix (ZMatrix.fits), + only used for the AEE model. delays : array of ints Array of MWA beamformer delay steps. Should be shape (n_pols, n_dipoles). Only applies to mwa_beam type files. @@ -4551,6 +4589,11 @@ def read( freq_range : array_like of float Range of frequencies to include in Hz, defaults to all available frequencies. Must be length 2. + mwa_include_cross_feed_coupling : bool + Option to include the couplings between the different feed directions + (i.e. couplings between x and y feeds). Including these couplings + is more correct, but the option is supplied to enable matching with + some codes that exclude these couplings. Default is True. Raises ------ @@ -4571,9 +4614,19 @@ def read( basename, extension = os.path.splitext(test_file) extension = extension.lower() if extension == ".fits" or extension == ".beamfits": - file_type = "beamfits" + from astropy.io import fits + + with fits.open(test_file) as fname: + if "CTYPE1" in fname[0].header: + file_type = "beamfits" + elif "FREQ" in fname[0].header: + file_type = "mwa_beam" + if mwa_model_type is None: + mwa_model_type = "aee" elif extension == ".hdf5" or extension == ".h5": file_type = "mwa_beam" + if mwa_model_type is None: + mwa_model_type = "fee" elif extension == ".txt" or extension == ".yaml": file_type = "cst" elif extension == ".ffe": @@ -4780,10 +4833,13 @@ def read( if file_type == "mwa_beam": self.read_mwa_beam( filename, + model_type=mwa_model_type, + zfile=mwa_zfile, delays=delays, amplitudes=amplitudes, pixels_per_deg=pixels_per_deg, freq_range=freq_range, + include_cross_feed_coupling=mwa_include_cross_feed_coupling, run_check=run_check, check_extra=check_extra, run_check_acceptability=run_check_acceptability, @@ -4964,22 +5020,40 @@ def _uvbeam_constructor(loader, node): if isinstance(values["filename"], str): files_use = [values["filename"]] + zfiles_use = None + if "mwa_zfile" in values: + zfiles_use = values["mwa_zfile"] + if isinstance(values["mwa_zfile"], str): + zfiles_use = [values["mwa_zfile"]] + if "path_variable" in values: path_var = values.pop("path_variable") # first check to see if this is a file on disk test_files = [os.path.join(path_var, file) for file in files_use] - files_exist = np.asarray([os.path.exists(file) for file in test_files]) + if "mwa_zfile" in values: + test_zfiles = [os.path.join(path_var, file) for file in zfiles_use] + else: + test_zfiles = [] + files_exist = np.asarray( + [os.path.exists(file) for file in (test_files + test_zfiles)] + ) if np.any(files_exist): files_use = test_files + zfiles_use = test_zfiles if not np.all(files_exist): - missing_files = (np.asarray(test_files))[np.nonzero(~files_exist)] - raise FileNotFoundError(f"File(s) {missing_files} do not exist.") + missing_files = (np.asarray(test_files + test_zfiles))[ + np.nonzero(~files_exist) + ] + raise FileNotFoundError( + f"File(s) {missing_files.tolist()} do not exist." + ) else: bad_pathvar_message = ( "If 'path_variable' is specified, it should either be the " "directory a file is in or take the form of a module.variable_name " "where the variable name can be imported from the module. " - f"The path_variable is {path_var}. Files {test_files} do not exist " + f"The path_variable is {path_var}. Files {test_files + test_zfiles} " + "do not exist " ) path_parts = (path_var).split(".") var_name = path_parts[-1] @@ -4999,11 +5073,21 @@ def _uvbeam_constructor(loader, node): ) from ie path_var = getattr(module, var_name) files_use = [os.path.join(path_var, file) for file in files_use] - files_exist = np.asarray([os.path.exists(file) for file in files_use]) + if "mwa_zfile" in values: + zfiles_use = [os.path.join(path_var, file) for file in zfiles_use] + else: + zfiles_use = [] + + files_exist = np.asarray( + [os.path.exists(file) for file in (files_use + zfiles_use)] + ) if not np.all(files_exist): - missing_files = (np.asarray(files_use))[np.nonzero(~files_exist)] + missing_files = (np.asarray(files_use + zfiles_use))[ + np.nonzero(~files_exist) + ] raise FileNotFoundError( - bad_pathvar_message + f"and file(s) {missing_files} do not exist." + bad_pathvar_message + f"and file(s) {missing_files.tolist()} " + "do not exist." ) for i, file in enumerate(files_use): @@ -5016,6 +5100,16 @@ def _uvbeam_constructor(loader, node): files_use = files_use[0] values["filename"] = files_use + if "mwa_zfile" in values: + for i, file in enumerate(zfiles_use): + # if file does not exist, check pyuvsim cache defined from astropy + # prescription treat file as download url to check astropy cache for file + if not os.path.exists(file) and is_url_in_cache(file, pkgname="pyuvsim"): + zfiles_use[i] = cache_contents("pyuvsim")[file] + if len(zfiles_use) == 1: + zfiles_use = zfiles_use[0] + values["mwa_zfile"] = zfiles_use + beam = UVBeam.from_file(**values) return beam diff --git a/tests/uvbeam/conftest.py b/tests/uvbeam/conftest.py index 501bb9f41a..e6b58d4a0d 100644 --- a/tests/uvbeam/conftest.py +++ b/tests/uvbeam/conftest.py @@ -334,7 +334,7 @@ def phased_array_beam_1freq(phased_array_beam_2freq): @pytest.fixture(scope="module") -def mwa_beam_1ppd_main(): +def mwa_fee_1ppd_main(): beam = UVBeam() beam.read_mwa_beam(fetch_data("mwa_full_EE"), pixels_per_deg=1) @@ -343,8 +343,51 @@ def mwa_beam_1ppd_main(): @pytest.fixture(scope="function") -def mwa_beam_1ppd(mwa_beam_1ppd_main): - beam = mwa_beam_1ppd_main.copy() +def mwa_fee_1ppd(mwa_fee_1ppd_main): + beam = mwa_fee_1ppd_main.copy() + + yield beam + del beam + + +@pytest.fixture(scope="module") +def mwa_aee_files(): + return {"jfile": fetch_data("mwa_jmatrix"), "zfile": fetch_data("mwa_zmatrix")} + + +@pytest.fixture(scope="module") +def mwa_aee_main(mwa_aee_files): + beam = UVBeam() + filename = mwa_aee_files["jfile"] + zfile = mwa_aee_files["zfile"] + beam.read_mwa_beam(filename, zfile=zfile) + + yield beam + del beam + + +@pytest.fixture(scope="function") +def mwa_aee(mwa_aee_main): + beam = mwa_aee_main.copy() + + yield beam + del beam + + +@pytest.fixture(scope="module") +def mwa_aee_noxy_main(mwa_aee_files): + beam = UVBeam() + filename = mwa_aee_files["jfile"] + zfile = mwa_aee_files["zfile"] + beam.read_mwa_beam(filename, zfile=zfile, include_cross_feed_coupling=False) + + yield beam + del beam + + +@pytest.fixture(scope="function") +def mwa_aee_noxy(mwa_aee_noxy_main): + beam = mwa_aee_noxy_main.copy() yield beam del beam diff --git a/tests/uvbeam/test_mwa_beam.py b/tests/uvbeam/test_mwa_beam.py index 1e3d46c806..f56f2e424c 100644 --- a/tests/uvbeam/test_mwa_beam.py +++ b/tests/uvbeam/test_mwa_beam.py @@ -1,8 +1,11 @@ # Copyright (c) 2019 Radio Astronomy Software Group # Licensed under the 2-clause BSD License +import copy +import re import numpy as np import pytest +from astropy.io import fits from pyuvdata import UVBeam, utils from pyuvdata.datasets import fetch_data @@ -10,27 +13,67 @@ from pyuvdata.uvbeam.mwa_beam import P1sin, P1sin_array -def test_read_write_mwa(mwa_beam_1ppd, tmp_path): +@pytest.fixture() +def mwabeam_kwargs(mwa_aee_files): + return { + "fee": {"pixels_per_deg": 1}, + "aee": {"zfile": mwa_aee_files["zfile"]}, + "aee_noxy": { + "zfile": mwa_aee_files["zfile"], + "include_cross_feed_coupling": False, + }, + } + + +@pytest.fixture() +def uvbeam_kwargs(mwabeam_kwargs): + uvbeam_kwargs = copy.deepcopy(mwabeam_kwargs) + uvbeam_kwargs["aee"]["mwa_zfile"] = uvbeam_kwargs["aee"].pop("zfile") + uvbeam_kwargs["aee_noxy"]["mwa_zfile"] = uvbeam_kwargs["aee_noxy"].pop("zfile") + uvbeam_kwargs["aee_noxy"]["mwa_include_cross_feed_coupling"] = uvbeam_kwargs[ + "aee_noxy" + ].pop("include_cross_feed_coupling") + return uvbeam_kwargs + + +@pytest.fixture() +def beam_filenames(mwa_aee_files): + return { + "fee": fetch_data("mwa_full_EE"), + "aee": mwa_aee_files["jfile"], + "aee_noxy": mwa_aee_files["jfile"], + } + + +@pytest.mark.parametrize("model", ["fee", "aee"]) +def test_read_write_mwa(beam_filenames, mwabeam_kwargs, model, tmp_path): """Basic read/write test.""" - beam1 = mwa_beam_1ppd - beam2 = UVBeam() - beam1.read_mwa_beam(fetch_data("mwa_full_EE"), pixels_per_deg=1) - assert beam1.filename == ["mwa_full_EE_test.h5"] + filename = beam_filenames[model] + kwargs = mwabeam_kwargs[model] + + beam1 = UVBeam() + beam2 = UVBeam() + beam1.read_mwa_beam(filename, **kwargs) + + if model == "fee": + assert beam1.filename == ["mwa_full_EE_test.h5"] + assert beam1.data_array.shape == (2, 2, 3, 91, 360) + # this is entirely empirical, just to prevent unexpected changes. + # The actual values have been validated through external tests against + # the mwa_pb repo. + assert np.isclose( + np.max(np.abs(beam1.data_array)), + 0.6823676193472403, + rtol=beam1._data_array.tols[0], + atol=beam1._data_array.tols[1], + ) + else: + assert beam1.filename == ["JMatrix_3freq.fits", "ZMatrix_3freq.fits"] + assert beam1.data_array.shape == (2, 2, 3, 31, 121) assert beam1.pixel_coordinate_system == "az_za" assert beam1.beam_type == "efield" - assert beam1.data_array.shape == (2, 2, 3, 91, 360) - - # this is entirely empirical, just to prevent unexpected changes. - # The actual values have been validated through external tests against - # the mwa_pb repo. - assert np.isclose( - np.max(np.abs(beam1.data_array)), - 0.6823676193472403, - rtol=beam1._data_array.tols[0], - atol=beam1._data_array.tols[1], - ) assert "x" in beam1.feed_array assert "y" in beam1.feed_array @@ -45,8 +88,22 @@ def test_read_write_mwa(mwa_beam_1ppd, tmp_path): @pytest.mark.filterwarnings("ignore:There are some terminated dipoles") -def test_mwa_orientation(mwa_beam_1ppd): - power_beam = mwa_beam_1ppd.efield_to_power(inplace=False) +@pytest.mark.parametrize("model", ["fee", "aee", "aee_noxy"]) +def test_mwa_orientation( + mwa_fee_1ppd, mwa_aee, mwa_aee_noxy, beam_filenames, uvbeam_kwargs, model +): + if model == "fee": + ebeam = mwa_fee_1ppd + near_hor_za = 80 + near_zen_za = 2 + else: + if model == "aee": + ebeam = mwa_aee + else: + ebeam = mwa_aee_noxy + near_hor_za = 81 + near_zen_za = 3 + power_beam = ebeam.efield_to_power(inplace=False) za_val = np.nonzero(np.isclose(power_beam.axis2_array, 15.0 * np.pi / 180)) @@ -58,13 +115,13 @@ def test_mwa_orientation(mwa_beam_1ppd): == utils.polstr2num( "ee", x_orientation=power_beam.get_x_orientation_from_feeds() ) - ) + )[0] north_ind = np.nonzero( power_beam.polarization_array == utils.polstr2num( "nn", x_orientation=power_beam.get_x_orientation_from_feeds() ) - ) + )[0] # check that the e/w dipole is more sensitive n/s assert ( @@ -85,11 +142,13 @@ def test_mwa_orientation(mwa_beam_1ppd): # single dipole delays = np.full((2, 16), 32, dtype=int) delays[:, 5] = 0 - efield_beam = UVBeam.from_file( - fetch_data("mwa_full_EE"), pixels_per_deg=1, delays=delays - ) - za_val = np.nonzero(np.isclose(efield_beam.axis2_array, 80.0 * np.pi / 180)) + filename = beam_filenames[model] + kwargs = uvbeam_kwargs[model] + kwargs["delays"] = delays + efield_beam = UVBeam.from_file(filename, **kwargs) + + za_val = np.nonzero(np.isclose(efield_beam.axis2_array, near_hor_za * np.pi / 180)) max_az_response = np.max(np.abs(efield_beam.data_array[0, east_ind, 0, za_val, :])) max_za_response = np.max(np.abs(efield_beam.data_array[1, east_ind, 0, za_val, :])) @@ -99,9 +158,10 @@ def test_mwa_orientation(mwa_beam_1ppd): max_za_response = np.max(np.abs(efield_beam.data_array[1, north_ind, 0, za_val, :])) assert max_az_response > max_za_response + # go back to zenith pointed full tile beam # check the sign of the responses are as expected close to zenith - efield_beam = mwa_beam_1ppd - za_val = np.nonzero(np.isclose(power_beam.axis2_array, 2.0 * np.pi / 180)) + efield_beam = ebeam + za_val = np.nonzero(np.isclose(power_beam.axis2_array, near_zen_za * np.pi / 180)) # first check zenith angle aligned response assert efield_beam.data_array[1, east_ind, 0, za_val, east_az] > 0 @@ -127,7 +187,10 @@ def test_mwa_orientation(mwa_beam_1ppd): ), ], ) -def test_mwa_pointing(delay_set, az_val, za_range): +@pytest.mark.parametrize("model", ["fee", "aee", "aee_noxy"]) +def test_mwa_pointing( + delay_set, az_val, za_range, beam_filenames, uvbeam_kwargs, model +): # Test that pointing the beam moves the peak in the right direction. delays = np.empty((2, 16), dtype=int) @@ -135,9 +198,10 @@ def test_mwa_pointing(delay_set, az_val, za_range): for pol in range(2): delays[pol] = delay_set - mwa_beam = UVBeam.from_file( - fetch_data("mwa_full_EE"), pixels_per_deg=1, beam_type="efield", delays=delays - ) + filename = beam_filenames[model] + kwargs = uvbeam_kwargs[model] + kwargs["delays"] = delays + mwa_beam = UVBeam.from_file(filename, **kwargs) mwa_beam.efield_to_power(calc_cross_pols=False) # set up zenith angle, azimuth and frequency arrays to evaluate with @@ -189,36 +253,44 @@ def test_mwa_pointing(delay_set, az_val, za_range): assert np.rad2deg(za_array[max_nn_loc]) < za_range[1] -def test_freq_range(mwa_beam_1ppd): - beam1 = mwa_beam_1ppd +@pytest.mark.parametrize("model", ["fee", "aee"]) +def test_freq_range(mwa_fee_1ppd, mwa_aee, model, beam_filenames, mwabeam_kwargs): + if model == "fee": + beam1 = mwa_fee_1ppd + else: + beam1 = mwa_aee + + filename = beam_filenames[model] + kwargs = mwabeam_kwargs[model] + beam2 = UVBeam() # include all - beam2.read_mwa_beam( - fetch_data("mwa_full_EE"), pixels_per_deg=1, freq_range=[100e6, 200e6] - ) + beam2.read_mwa_beam(filename, **kwargs, freq_range=[100e6, 200e6]) assert beam1 == beam2 - beam2.read_mwa_beam( - fetch_data("mwa_full_EE"), pixels_per_deg=1, freq_range=[100e6, 150e6] - ) + beam2.read_mwa_beam(filename, **kwargs, freq_range=[100e6, 170e6]) beam1.select(freq_chans=[0, 1]) assert beam1.history != beam2.history beam1.history = beam2.history assert beam1 == beam2 + +def test_freq_range_errors(): + beam1 = UVBeam() + with check_warnings(UserWarning, match="Only one available frequency"): beam1.read_mwa_beam( fetch_data("mwa_full_EE"), pixels_per_deg=1, freq_range=[100e6, 130e6] ) with pytest.raises(ValueError, match="No frequencies available in freq_range"): - beam2.read_mwa_beam( + beam1.read_mwa_beam( fetch_data("mwa_full_EE"), pixels_per_deg=1, freq_range=[100e6, 110e6] ) with pytest.raises(ValueError, match="freq_range must have 2 elements."): - beam2.read_mwa_beam( + beam1.read_mwa_beam( fetch_data("mwa_full_EE"), pixels_per_deg=1, freq_range=[100e6] ) @@ -245,31 +317,27 @@ def test_bad_amps(): beam1 = UVBeam() amps = np.ones([2, 8]) - with pytest.raises(ValueError) as cm: + with pytest.raises(ValueError, match="amplitudes must be shape"): beam1.read_mwa_beam( fetch_data("mwa_full_EE"), pixels_per_deg=1, amplitudes=amps ) - assert str(cm.value).startswith("amplitudes must be shape") def test_bad_delays(): beam1 = UVBeam() delays = np.zeros([2, 8], dtype="int") - with pytest.raises(ValueError) as cm: + with pytest.raises(ValueError, match="delays must be shape"): beam1.read_mwa_beam(fetch_data("mwa_full_EE"), pixels_per_deg=1, delays=delays) - assert str(cm.value).startswith("delays must be shape") delays = np.zeros((2, 16), dtype="int") delays = delays + 64 - with pytest.raises(ValueError) as cm: + with pytest.raises(ValueError, match="There are delays greater than 32"): beam1.read_mwa_beam(fetch_data("mwa_full_EE"), pixels_per_deg=1, delays=delays) - assert str(cm.value).startswith("There are delays greater than 32") delays = np.zeros((2, 16), dtype="float") - with pytest.raises(ValueError) as cm: + with pytest.raises(ValueError, match="Delays must be integers."): beam1.read_mwa_beam(fetch_data("mwa_full_EE"), pixels_per_deg=1, delays=delays) - assert str(cm.value).startswith("Delays must be integers.") def test_dead_dipoles(): @@ -301,3 +369,97 @@ def test_dead_dipoles(): + beam1.pyuvdata_version_str ) assert utils.history._check_histories(history_str, beam1.history) + + +@pytest.mark.parametrize( + ("change", "errmsg"), + [ + ( + "diff_exten", + "model_type could not be determined for MWA beam file, use " + "the model_type keyword to specify the type (one of 'fee' " + "or 'aee'). Filename is", + ), + ("no_grid", "Data does not appear to be on a grid"), + ("bad_th_reshape", "thetas do not appear to be on expected grid"), + ("bad_ph_reshape", "phis do not appear to be on expected grid"), + ("diff_th", "Inconsistent theta values across frequencies"), + ("diff_ph", "Inconsistent phi values across frequencies"), + ], +) +def test_aee_jfile_errors(tmp_path, mwa_aee_files, change, errmsg): + + if change == "diff_exten": + testfile = str(tmp_path / "test_aee_j.foo") + else: + testfile = str(tmp_path / "test_aee_j.fits") + + with fits.open(mwa_aee_files["jfile"]) as jfile: + first_hdu = jfile[0] + data = first_hdu.data + if change == "no_grid": + # change the first theta value slightly + data[0, 0] += 0.1 + elif change == "bad_th_reshape": + # set the first theta to 2nd theta + data[0, 0] = data[1, 0] + elif change == "bad_ph_reshape": + # set the first phi to 2nd phi + data[0, 1] = data[100, 1] + elif change == "diff_th": + # set the first phi to 2nd phi + data[:, 0] += 1 + elif change == "diff_ph": + # set the first phi to 2nd phi + data[:, 1] += 1 + + first_hdu.data = data + hdulist = [first_hdu] + for nhdu in range(1, len(jfile)): + hdulist.append(jfile[nhdu]) + + hdulist = fits.HDUList(hdulist) + + hdulist.writeto(testfile, overwrite=True) + hdulist.close() + + with pytest.raises(ValueError, match=re.escape(errmsg)): + UVBeam.from_file( + testfile, mwa_zfile=mwa_aee_files["zfile"], file_type="mwa_beam" + ) + + +@pytest.mark.parametrize( + ("change", "errmsg"), + [ + ( + "diff_nf", + "Zmatrix file does not have as the same number of frequencies as " + "Jmatrix file.", + ), + ("diff_f", "Zmatrix 0th freq does not match Jmatrix file."), + ], +) +def test_aee_zfile_errors(tmp_path, mwa_aee_files, change, errmsg): + + testfile = str(tmp_path / "test_aee_z.fits") + + with fits.open(mwa_aee_files["zfile"]) as zfile: + first_hdu = zfile[0] + n_freqs = len(zfile) + if change == "diff_nf": + n_freqs = n_freqs - 1 + elif change == "diff_f": + first_hdu.header["freq"] += 2 + + hdulist = [first_hdu] + for nhdu in range(1, n_freqs): + hdulist.append(zfile[nhdu]) + + hdulist = fits.HDUList(hdulist) + + hdulist.writeto(testfile, overwrite=True) + hdulist.close() + + with pytest.raises(ValueError, match=errmsg): + UVBeam.from_file(mwa_aee_files["jfile"], mwa_zfile=testfile) diff --git a/tests/uvbeam/test_uvbeam.py b/tests/uvbeam/test_uvbeam.py index b2228f321b..d9e110009a 100644 --- a/tests/uvbeam/test_uvbeam.py +++ b/tests/uvbeam/test_uvbeam.py @@ -3026,17 +3026,22 @@ def test_generic_read_cst(): @pytest.mark.filterwarnings("ignore:This beamfits file has no information about") @pytest.mark.filterwarnings("ignore:The mount_type keyword is set") @pytest.mark.parametrize( - "dname", ["hera_fagnoni_dipole_yaml", "mwa_full_EE", "hera_casa_beam"] + "dname", ["hera_fagnoni_dipole_yaml", "mwa_full_EE", "mwa_AEE", "hera_casa_beam"] ) -def test_generic_read(dname): +def test_generic_read(mwa_aee_files, dname): """Test generic read can infer the file types correctly.""" uvb = UVBeam() # going to check in a second anyway, no need to double check. - filename = fetch_data(dname) + if dname == "mwa_AEE": + filename = mwa_aee_files["jfile"] + zfile = mwa_aee_files["zfile"] + else: + filename = fetch_data(dname) + zfile = None if "yaml" in dname: # make sure the datafiles are cached fetch_data(["hera_fagnoni_dipole_150", "hera_fagnoni_dipole_123"]) - uvb.read(filename, mount_type="fixed", run_check=False) + uvb.read(filename, mount_type="fixed", mwa_zfile=zfile, run_check=False) # hera casa beam is missing some parameters but we just want to check # that reading is going okay if dname == "hera_casa_beam": @@ -3151,18 +3156,25 @@ def test_generic_read_all_bad_files(tmp_path): @pytest.mark.filterwarnings("ignore:This beamfits file has no information about") @pytest.mark.filterwarnings("ignore:The mount_type keyword is set") @pytest.mark.parametrize( - "dname", ["hera_fagnoni_dipole_yaml", "mwa_full_EE", "hera_casa_beam"] + "dname", ["hera_fagnoni_dipole_yaml", "mwa_full_EE", "mwa_AEE", "hera_casa_beam"] ) -def test_from_file(dname): +def test_from_file(mwa_aee_files, dname): """Test from file produces same the results as reading explicitly.""" uvb = UVBeam() # don't run checks because of casa_beamfits, we'll do that later - filename = fetch_data(dname) + if dname == "mwa_AEE": + filename = mwa_aee_files["jfile"] + zfile = mwa_aee_files["zfile"] + else: + filename = fetch_data(dname) + zfile = None if "yaml" in dname: # make sure the datafiles are cached fetch_data(["hera_fagnoni_dipole_150", "hera_fagnoni_dipole_123"]) - uvb2 = UVBeam.from_file(filename, mount_type="fixed", run_check=False) - uvb.read(filename, mount_type="fixed", run_check=False) + uvb2 = UVBeam.from_file( + filename, mount_type="fixed", mwa_zfile=zfile, run_check=False + ) + uvb.read(filename, mount_type="fixed", mwa_zfile=zfile, run_check=False) # hera casa beam is missing some parameters but we just want to check # that reading is going okay if dname == "hera_casa_beam": @@ -3192,11 +3204,22 @@ def test_from_file(dname): [ ["hera_fagnoni_dipole_yaml", True, False, True], ["mwa_full_EE", True, True, False], + ["mwa_AEE", True, True, False], + ["mwa_AEE", False, False, False], ["hera_casa_beam", False, False, True], ], ) -def test_yaml_constructor(dname, path_var, relative_path, file_list, tmp_path): - filename = fetch_data(dname) +def test_yaml_constructor( + mwa_aee_files, dname, path_var, relative_path, file_list, tmp_path +): + if dname == "mwa_AEE": + filename = mwa_aee_files["jfile"] + zfile = mwa_aee_files["zfile"] + include_cross_feed_coupling = path_var + else: + filename = fetch_data(dname) + zfile = None + include_cross_feed_coupling = True if "yaml" in dname: # make sure the datafiles are cached fetch_data(["hera_fagnoni_dipole_150", "hera_fagnoni_dipole_123"]) @@ -3223,11 +3246,26 @@ def test_yaml_constructor(dname, path_var, relative_path, file_list, tmp_path): run_check: False freq_range: [110.0e+6, 190.0e+6] """ + if zfile is not None: + if path_var: + dirname, zfile_use = os.path.split(zfile) + else: + zfile_use = zfile + input_yaml += f""" + mwa_zfile: {zfile_use} + mwa_include_cross_feed_coupling: {include_cross_feed_coupling} + """ beam_from_yaml = yaml.safe_load(input_yaml)["beam"] # don't run checks because of casa_beamfits, we'll do that later - uvb = UVBeam.from_file(filename, mount_type="fixed", run_check=False) + uvb = UVBeam.from_file( + filename, + mount_type="fixed", + mwa_zfile=zfile, + mwa_include_cross_feed_coupling=include_cross_feed_coupling, + run_check=False, + ) # hera casa beam is missing some parameters but we just want to check # that reading is going okay if dname == "hera_casa_beam": @@ -3276,15 +3314,22 @@ def test_yaml_constructor(dname, path_var, relative_path, file_list, tmp_path): output_yaml = yaml.safe_dump({"beam": uvb}) -def test_yaml_constructor_for_cached_file(): +@pytest.mark.parametrize("model", ["fee", "aee"]) +def test_yaml_constructor_for_cached_file(mwa_aee_files, model): """test that the yaml constructor uses pyuvsim astropy cache as a fallback""" - # url / key to access the cached file - dummy_url = "file://mwa_beam_file.h5" + if model == "fee": + filename = fetch_data("mwa_full_EE") + # url / key to access the cached file + dummy_url = "file://mwa_beam_file.h5" + else: + filename = mwa_aee_files["jfile"] + zfile = mwa_aee_files["zfile"] + # url / key to access the cached file + dummy_url = "file://JMatrix.fits" + dummy_url_zfile = "file://ZMatrix.fits" # load mwa uvbeam file from the data path to the astropy cache for pyuvsim - import_file_to_cache( - dummy_url, fetch_data("mwa_full_EE"), remove_original=False, pkgname="pyuvsim" - ) + import_file_to_cache(dummy_url, filename, remove_original=False, pkgname="pyuvsim") # check that file has been imported to cache assert is_url_in_cache(dummy_url, pkgname="pyuvsim") @@ -3292,43 +3337,87 @@ def test_yaml_constructor_for_cached_file(): # get full filepath for comparison filepath = cache_contents("pyuvsim")[dummy_url] + if model == "aee": + import_file_to_cache( + dummy_url_zfile, zfile, remove_original=False, pkgname="pyuvsim" + ) + assert is_url_in_cache(dummy_url_zfile, pkgname="pyuvsim") + zfilepath = cache_contents("pyuvsim")[dummy_url_zfile] + else: + zfilepath = None + # create yaml to safeload - input_yaml = f""" - beam: !UVBeam - filename: {dummy_url} - file_type: mwa_beam - """ + if model == "fee": + input_yaml = f""" + beam: !UVBeam + filename: {dummy_url} + file_type: mwa_beam + mwa_model_type: {model} + """ + else: + input_yaml = f""" + beam: !UVBeam + filename: {dummy_url} + file_type: mwa_beam + mwa_model_type: {model} + mwa_zfile: {dummy_url_zfile} + """ # compare loading uvbeam directly from cache and using safe_load to confirm # uvbeam loads properly # safe_load the beam beam_from_yaml = yaml.safe_load(input_yaml)["beam"] # load the beam with from_file - beam_from_file = UVBeam.from_file(filepath, file_type="mwa_beam") + beam_from_file = UVBeam.from_file( + filepath, file_type="mwa_beam", mwa_model_type=model, mwa_zfile=zfilepath + ) assert beam_from_yaml == beam_from_file -def test_yaml_constructor_errors(): - fname_use = fetch_data("hera_fagnoni_dipole_yaml") - dirname, basename = os.path.split(fname_use) +def test_yaml_constructor_errors(mwa_aee_files): + hera_fname_use = fetch_data("hera_fagnoni_dipole_yaml") + hera_dirname, hera_basename = os.path.split(hera_fname_use) + + mwaj_fname_use = mwa_aee_files["jfile"] + mwaj_dirname, mwaj_basename = os.path.split(mwaj_fname_use) + mwaz_fname_use = mwa_aee_files["zfile"] + mwaz_dirname, mwaz_basename = os.path.split(mwaz_fname_use) input_yaml = f""" beam: !UVBeam - filename: [{basename}, foo.yaml] - path_variable: {dirname} + filename: [{hera_basename}, foo.yaml] + path_variable: {hera_dirname} run_check: False """ with pytest.raises( FileNotFoundError, - match=re.escape(f"File(s) {[os.path.join(dirname, 'foo.yaml')]} do not exist."), + match=re.escape( + f"File(s) {[os.path.join(hera_dirname, 'foo.yaml')]} do not exist." + ), ): yaml.safe_load(input_yaml)["beam"] input_yaml = f""" beam: !UVBeam - filename: {basename} + filename: [{mwaj_basename}] + path_variable: {mwaj_dirname} + mwa_zfile: {"foo.fits"} + run_check: False + """ + + with pytest.raises( + FileNotFoundError, + match=re.escape( + f"File(s) {[os.path.join(mwaj_dirname, 'foo.fits')]} do not exist." + ), + ): + yaml.safe_load(input_yaml)["beam"] + + input_yaml = f""" + beam: !UVBeam + filename: {hera_basename} path_variable: DATA_PATH run_check: False """ @@ -3339,34 +3428,62 @@ def test_yaml_constructor_errors(): "If 'path_variable' is specified, it should either be the directory " "a file is in or take the form of a module.variable_name where the " "variable name can be imported from the module. The path_variable is " - f"DATA_PATH. Files {[os.path.join('DATA_PATH', basename)]} do not exist " - "and the path_variable does not have the form of a module.variable_name." + f"DATA_PATH. Files {[os.path.join('DATA_PATH', hera_basename)]} do " + "not exist and the path_variable does not have the form of a " + "module.variable_name." + ), + ): + yaml.safe_load(input_yaml)["beam"] + + input_yaml = f""" + beam: !UVBeam + filename: {hera_basename} + path_variable: pyuvdata.data.DATA_PATH + run_check: False + """ + bad_fname1 = [os.path.join("pyuvdata.data.DATA_PATH", hera_basename)] + bad_fname2 = [os.path.join(DATA_PATH, hera_basename)] + with pytest.raises( + FileNotFoundError, + match=re.escape( + "If 'path_variable' is specified, it should either be the directory " + "a file is in or take the form of a module.variable_name where the " + "variable name can be imported from the module. The path_variable is " + f"pyuvdata.data.DATA_PATH. Files {bad_fname1} do not exist and " + f"file(s) {bad_fname2} do not exist." ), ): yaml.safe_load(input_yaml)["beam"] input_yaml = f""" beam: !UVBeam - filename: {basename} + filename: {mwaj_basename} + mwa_zfile: {mwaz_basename} path_variable: pyuvdata.data.DATA_PATH run_check: False """ + bad_fnames1 = [ + os.path.join("pyuvdata.data.DATA_PATH", file) + for file in [mwaj_basename, mwaz_basename] + ] + bad_fnames2 = [ + os.path.join(DATA_PATH, file) for file in [mwaj_basename, mwaz_basename] + ] with pytest.raises( FileNotFoundError, match=re.escape( "If 'path_variable' is specified, it should either be the directory " "a file is in or take the form of a module.variable_name where the " "variable name can be imported from the module. The path_variable is " - "pyuvdata.data.DATA_PATH. Files " - f"{[os.path.join('pyuvdata.data.DATA_PATH', basename)]} do not exist " - f"and file(s) {[os.path.join(DATA_PATH, 'NicCSTbeams.yaml')]} do not exist." + f"pyuvdata.data.DATA_PATH. Files {bad_fnames1} do not exist and " + f"file(s) {bad_fnames2} do not exist." ), ): yaml.safe_load(input_yaml)["beam"] input_yaml = f""" beam: !UVBeam - filename: {basename} + filename: {hera_basename} path_variable: foo.DATA_PATH run_check: False """ @@ -3377,7 +3494,7 @@ def test_yaml_constructor_errors(): "a file is in or take the form of a module.variable_name where the " "variable name can be imported from the module. The path_variable is " "foo.DATA_PATH. Files " - f"{[os.path.join('foo.DATA_PATH', basename)]} do not exist " + f"{[os.path.join('foo.DATA_PATH', hera_basename)]} do not exist " f"and the module foo is not importable." ), ): @@ -3396,7 +3513,7 @@ def test_yaml_constructor_errors(): input_yaml = f""" beam: !UVBeam - filename: {fetch_data("hera_fagnoni_dipole_yaml")} + filename: {hera_fname_use} mount_type: fixed run_check: False freq_range: [120.0e+6] @@ -3516,7 +3633,7 @@ def test_fix_feeds_dep_warnings(cst_power_2freq_cut, mod_params, warn_msg): def test_plotting( tmp_path, cst_efield_1freq, - mwa_beam_1ppd, + mwa_fee_1ppd, telescope, uvbfuncs, uvbfuncs_kwargs, @@ -3535,7 +3652,7 @@ def test_plotting( if telescope == "hera": beam = cst_efield_1freq elif telescope == "mwa": - beam = mwa_beam_1ppd + beam = mwa_fee_1ppd if uvbfuncs is not None: if "to_healpix" in uvbfuncs: @@ -3554,8 +3671,8 @@ def test_plotting( ) -def test_plotting_errors(mwa_beam_1ppd): - beam = mwa_beam_1ppd +def test_plotting_errors(mwa_fee_1ppd): + beam = mwa_fee_1ppd if importlib.util.find_spec("matplotlib") is None: with pytest.raises(