Skip to content

Commit

Permalink
Merge pull request #174 from AstarVienna/fh/patchnb
Browse files Browse the repository at this point in the history 
Patch problematic notebook
  • Loading branch information
@teutoburg
teutoburg committed Apr 19, 2024
2 parents fd52794 + 96d5372 commit 4dd3b1b
Showing 1 changed file with 116 additions and 109 deletions.
225 changes: 116 additions & 109 deletions METIS/docs/example_notebooks/LSS-YSO_model_simulation.ipynb
View file
Expand Up @@ -282,13 +282,13 @@
"metadata": {},
"outputs": [],
"source": [
"meta = metis['spectral_traces'].meta\n",
"# meta = metis['spectral_traces'].meta\n",
"\n",
"det_xi = meta['CRVAL1'] + meta['CDELT1'] * (np.arange(2048) + 1 - meta['CRPIX1']) * u.Unit(meta['CUNIT1'])\n",
"det_xi = det_xi.to(u.arcsec)\n",
"# det_xi = meta['CRVAL1'] + meta['CDELT1'] * (np.arange(2048) + 1 - meta['CRPIX1']) * u.Unit(meta['CUNIT1'])\n",
"# det_xi = det_xi.to(u.arcsec)\n",
"\n",
"det_wave = (meta['CRVAL2'] + meta['CDELT2'] * (np.arange(2048) + 1 - meta['CRPIX2'])) * u.Unit(meta['CUNIT2'])\n",
"det_wave = det_wave.to(u.um) "
"# det_wave = (meta['CRVAL2'] + meta['CDELT2'] * (np.arange(2048) + 1 - meta['CRPIX2'])) * u.Unit(meta['CUNIT2'])\n",
"# det_wave = det_wave.to(u.um) "
]
},
{
Expand All @@ -297,10 +297,10 @@
"metadata": {},
"outputs": [],
"source": [
"ilim = np.array([750, 1300]) # pixels along spatial direction\n",
"jlim = np.array([300, 1750]) # pixels along wavelength direction\n",
"xlim = det_xi[ilim].value\n",
"ylim = det_wave[jlim].value"
"# ilim = np.array([750, 1300]) # pixels along spatial direction\n",
"# jlim = np.array([300, 1750]) # pixels along wavelength direction\n",
"# xlim = det_xi[ilim].value\n",
"# ylim = det_wave[jlim].value"
]
},
{
Expand All @@ -309,14 +309,14 @@
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(15, 8))\n",
"for i, (model, result) in enumerate(results.items()):\n",
" plt.subplot(1, 3, i+1)\n",
" plt.imshow((result[1].data - bgresult[1].data)[jlim[0]:jlim[1], ilim[0]:ilim[1]], origin='lower', norm=LogNorm(),\n",
" extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
" plt.xlabel(r\"position along slit [arcsec]\")\n",
" plt.ylabel(r\"Wavelength [$\\mu$m]\")\n",
" plt.title(model);"
"# plt.figure(figsize=(15, 8))\n",
"# for i, (model, result) in enumerate(results.items()):\n",
"# plt.subplot(1, 3, i+1)\n",
"# plt.imshow((result[1].data - bgresult[1].data)[jlim[0]:jlim[1], ilim[0]:ilim[1]], origin='lower', norm=LogNorm(),\n",
"# extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
"# plt.xlabel(r\"position along slit [arcsec]\")\n",
"# plt.ylabel(r\"Wavelength [$\\mu$m]\")\n",
"# plt.title(model);"
]
},
{
Expand All @@ -332,12 +332,12 @@
"metadata": {},
"outputs": [],
"source": [
"from pathlib import Path\n",
"for i, (model, result) in enumerate(results.items()):\n",
" outfile = Path(fitsfiles[model]).stem + \"-simulated_LSS_L\" + Path(fitsfiles[model]).suffix\n",
" result.writeto(outfile, overwrite=True)\n",
" print(fitsfiles[model], \"--->\", outfile)\n",
"bgresult.writeto(\"models_Lband_HD100546-background_simulated_LSS_L.fits\", overwrite=True)"
"# from pathlib import Path\n",
"# for i, (model, result) in enumerate(results.items()):\n",
"# outfile = Path(fitsfiles[model]).stem + \"-simulated_LSS_L\" + Path(fitsfiles[model]).suffix\n",
"# result.writeto(outfile, overwrite=True)\n",
"# print(fitsfiles[model], \"--->\", outfile)\n",
"# bgresult.writeto(\"models_Lband_HD100546-background_simulated_LSS_L.fits\", overwrite=True)"
]
},
{
Expand Down Expand Up @@ -457,15 +457,15 @@
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(15, 8))\n",
"for i, (model, result) in enumerate(rotresults.items()):\n",
" plt.subplot(1, 3, i+1)\n",
" plt.imshow((result[1].data - bgresult[1].data)[jlim[0]:jlim[1], ilim[0]:ilim[1]], \n",
" origin='lower', norm=LogNorm(),\n",
" extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
" plt.xlabel(r\"position along slit [arcsec]\")\n",
" plt.ylabel(r\"Wavelength [$\\mu$m]\")\n",
" plt.title(model + \", rotated by \" + str(angle) + \" degrees\");"
"# plt.figure(figsize=(15, 8))\n",
"# for i, (model, result) in enumerate(rotresults.items()):\n",
"# plt.subplot(1, 3, i+1)\n",
"# plt.imshow((result[1].data - bgresult[1].data)[jlim[0]:jlim[1], ilim[0]:ilim[1]], \n",
"# origin='lower', norm=LogNorm(),\n",
"# extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
"# plt.xlabel(r\"position along slit [arcsec]\")\n",
"# plt.ylabel(r\"Wavelength [$\\mu$m]\")\n",
"# plt.title(model + \", rotated by \" + str(angle) + \" degrees\")"
]
},
{
Expand All @@ -474,10 +474,10 @@
"metadata": {},
"outputs": [],
"source": [
"for i, (model, result) in enumerate(rotresults.items()):\n",
" outfile = Path(fitsfiles[model]).stem + \"-rot_\" + str(angle) + \"-simulated_LSS_L\" + Path(fitsfiles[model]).suffix\n",
" result.writeto(outfile, overwrite=True)\n",
" print(fitsfiles[model], \"--->\", outfile)"
"# for i, (model, result) in enumerate(rotresults.items()):\n",
"# outfile = Path(fitsfiles[model]).stem + \"-rot_\" + str(angle) + \"-simulated_LSS_L\" + Path(fitsfiles[model]).suffix\n",
"# result.writeto(outfile, overwrite=True)\n",
"# print(fitsfiles[model], \"--->\", outfile)"
]
},
{
Expand All @@ -493,21 +493,21 @@
"metadata": {},
"outputs": [],
"source": [
"fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(6, 12))\n",
"# fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1, figsize=(6, 12))\n",
"\n",
"ax1.plot(det_xi, (results['emptycav'][1].data - bgresult[1].data)[500, :], label=\"no rotation\")\n",
"ax1.plot(det_xi, (rotresults['emptycav'][1].data - bgresult[1].data)[500, :], label=f\"rotation {angle} deg\")\n",
"ax1.set_xlim(xlim[0], xlim[-1])\n",
"ax1.set_ylim(2e4, 1e9)\n",
"ax1.set_xlabel(\"arcsec\")\n",
"ax1.semilogy()\n",
"ax1.legend();\n",
"# ax1.plot(det_xi, (results['emptycav'][1].data - bgresult[1].data)[500, :], label=\"no rotation\")\n",
"# ax1.plot(det_xi, (rotresults['emptycav'][1].data - bgresult[1].data)[500, :], label=f\"rotation {angle} deg\")\n",
"# ax1.set_xlim(xlim[0], xlim[-1])\n",
"# ax1.set_ylim(2e4, 1e9)\n",
"# ax1.set_xlabel(\"arcsec\")\n",
"# ax1.semilogy()\n",
"# ax1.legend()\n",
"\n",
"fig.subplots_adjust(left=0.1)\n",
"ax2.imshow(rotsources['emptycav'].cube_fields[0].data.sum(axis=0) + 1e-14, norm=LogNorm(vmin=1e-4, vmax=10),\n",
" extent=(x_lim[0], x_lim[-1], y_lim[0], y_lim[-1]))\n",
"ax2.set_xlabel(\"arcsec\")\n",
"ax2.set_ylabel(\"arcsec\");"
"# fig.subplots_adjust(left=0.1)\n",
"# ax2.imshow(rotsources['emptycav'].cube_fields[0].data.sum(axis=0) + 1e-14, norm=LogNorm(vmin=1e-4, vmax=10),\n",
"# extent=(x_lim[0], x_lim[-1], y_lim[0], y_lim[-1]))\n",
"# ax2.set_xlabel(\"arcsec\")\n",
"# ax2.set_ylabel(\"arcsec\")"
]
},
{
Expand Down Expand Up @@ -573,21 +573,21 @@
"metadata": {},
"outputs": [],
"source": [
"xmin, xmax = 975, 1075\n",
"xaxis = np.arange(xmin, xmax)\n",
"y_1, y_2 = 200, 1800\n",
"lam_1, lam_2 = det_wave[y_1], det_wave[y_2]\n",
"# xmin, xmax = 975, 1075\n",
"# xaxis = np.arange(xmin, xmax)\n",
"# y_1, y_2 = 200, 1800\n",
"# lam_1, lam_2 = det_wave[y_1], det_wave[y_2]\n",
"\n",
"plt.figure(figsize=(10, 5))\n",
"plt.plot(xaxis, std_bgsub[200, xmin:xmax], label=fr\"$\\lambda = {lam_1:.2f}\\,\\mu\\mathrm{{m}}$\")\n",
"plt.plot(xaxis, std_bgsub[1800, xmin:xmax], label=fr\"$\\lambda = {lam_2:.2f}\\,\\mu\\mathrm{{m}}$\")\n",
"plt.plot(xaxis, std_bgsub[:, xmin:xmax].mean(axis=0), label='average')\n",
"plt.vlines((1024 - 10, 1024 + 10), 1, 1e5, colors='black', linestyles='dashed', label='extraction aperture')\n",
"plt.xlabel(\"pixel\")\n",
"plt.ylabel(\"e/s\")\n",
"plt.semilogy()\n",
"plt.legend()\n",
"plt.title(\"spatial cuts through standard spectrum\");"
"# plt.figure(figsize=(10, 5))\n",
"# plt.plot(xaxis, std_bgsub[200, xmin:xmax], label=fr\"$\\lambda = {lam_1:.2f}\\,\\mu\\mathrm{{m}}$\")\n",
"# plt.plot(xaxis, std_bgsub[1800, xmin:xmax], label=fr\"$\\lambda = {lam_2:.2f}\\,\\mu\\mathrm{{m}}$\")\n",
"# plt.plot(xaxis, std_bgsub[:, xmin:xmax].mean(axis=0), label='average')\n",
"# plt.vlines((1024 - 10, 1024 + 10), 1, 1e5, colors='black', linestyles='dashed', label='extraction aperture')\n",
"# plt.xlabel(\"pixel\")\n",
"# plt.ylabel(\"e/s\")\n",
"# plt.semilogy()\n",
"# plt.legend()\n",
"# plt.title(\"spatial cuts through standard spectrum\");"
]
},
{
Expand Down Expand Up @@ -633,11 +633,11 @@
"metadata": {},
"outputs": [],
"source": [
"plt.plot(det_wave, std_1d)\n",
"plt.xlabel('Wavelength (um)')\n",
"exptime = std_result[1].header['EXPTIME']\n",
"plt.ylabel('electrons/s')\n",
"plt.title(r\"Extracted standard spectrum ($T_{\\mathrm{exp}} = 1\\,\\mathrm{s}$)\");"
"# plt.plot(det_wave, std_1d)\n",
"# plt.xlabel('Wavelength (um)')\n",
"# exptime = std_result[1].header['EXPTIME']\n",
"# plt.ylabel('electrons/s')\n",
"# plt.title(r\"Extracted standard spectrum ($T_{\\mathrm{exp}} = 1\\,\\mathrm{s}$)\");"
]
},
{
Expand All @@ -653,12 +653,12 @@
"metadata": {},
"outputs": [],
"source": [
"flux_conv = 1 * u.Jy / (std_1d * u.electron)\n",
"plt.plot(det_wave, flux_conv)\n",
"plt.ylim(0, 1e-5)\n",
"plt.xlabel(\"Wavelength (um)\")\n",
"plt.ylabel(\"Jy / (e/s)\")\n",
"plt.title(\"Flux conversion\");"
"# flux_conv = 1 * u.Jy / (std_1d * u.electron)\n",
"# plt.plot(det_wave, flux_conv)\n",
"# plt.ylim(0, 1e-5)\n",
"# plt.xlabel(\"Wavelength (um)\")\n",
"# plt.ylabel(\"Jy / (e/s)\")\n",
"# plt.title(\"Flux conversion\");"
]
},
{
Expand All @@ -674,8 +674,8 @@
"metadata": {},
"outputs": [],
"source": [
"emptycav_Jy = flux_conv[:, None] * (results['emptycav'][1].data - bgresult[1].data) / exptime\n",
"emptycav_Jy_rot = flux_conv[:, None] * (rotresults['emptycav'][1].data - bgresult[1].data) / exptime"
"# emptycav_Jy = flux_conv[:, None] * (results['emptycav'][1].data - bgresult[1].data) / exptime\n",
"# emptycav_Jy_rot = flux_conv[:, None] * (rotresults['emptycav'][1].data - bgresult[1].data) / exptime"
]
},
{
Expand All @@ -684,22 +684,22 @@
"metadata": {},
"outputs": [],
"source": [
"fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(13, 12))\n",
"ax1.imshow(emptycav_Jy[jlim[0]:jlim[1], ilim[0]:ilim[1]].value, origin='lower', \n",
" norm=LogNorm(vmin=1e-5, vmax=1), extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
"ax1.set_xlabel(r\"position along slit [arcsec]\")\n",
"ax1.set_ylabel(r\"Wavelength [$\\mu$m]\") \n",
"ax1.set_title(model + \", no rotation\")\n",
"# fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(13, 12))\n",
"# ax1.imshow(emptycav_Jy[jlim[0]:jlim[1], ilim[0]:ilim[1]].value, origin='lower', \n",
"# norm=LogNorm(vmin=1e-5, vmax=1), extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
"# ax1.set_xlabel(r\"position along slit [arcsec]\")\n",
"# ax1.set_ylabel(r\"Wavelength [$\\mu$m]\") \n",
"# ax1.set_title(model + \", no rotation\")\n",
"\n",
"img = ax2.imshow(emptycav_Jy_rot[jlim[0]:jlim[1], ilim[0]:ilim[1]].value, origin='lower',\n",
" norm=LogNorm(vmin=1e-5, vmax=1), extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
"ax2.set_xlabel(r\"position along slit [arcsec]\")\n",
"ax2.set_ylabel(r\"Wavelength [$\\mu$m]\") \n",
"ax2.set_title(model + \", rotated by \" + str(angle) + \" degrees\");\n",
"# img = ax2.imshow(emptycav_Jy_rot[jlim[0]:jlim[1], ilim[0]:ilim[1]].value, origin='lower',\n",
"# norm=LogNorm(vmin=1e-5, vmax=1), extent=(xlim[0], xlim[1], ylim[0], ylim[1]), aspect='auto')\n",
"# ax2.set_xlabel(r\"position along slit [arcsec]\")\n",
"# ax2.set_ylabel(r\"Wavelength [$\\mu$m]\") \n",
"# ax2.set_title(model + \", rotated by \" + str(angle) + \" degrees\");\n",
"\n",
"fig.subplots_adjust(right=0.8)\n",
"cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])\n",
"fig.colorbar(img, cax=cbar_ax, label='Jy')"
"# fig.subplots_adjust(right=0.8)\n",
"# cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])\n",
"# fig.colorbar(img, cax=cbar_ax, label='Jy')"
]
},
{
Expand All @@ -715,21 +715,21 @@
"metadata": {},
"outputs": [],
"source": [
"for i, (model, result) in enumerate(results.items()):\n",
" result[1].data = flux_conv[:, None].value * (result[1].data - bgresult[1].data) / exptime\n",
" result[1].header['BUNIT'] = \"e/s\"\n",
" outfile = (Path(fitsfiles[model]).stem + \"-simulated_LSS_L-bgsub_fluxcal\" \n",
" + Path(fitsfiles[model]).suffix)\n",
" result.writeto(outfile, overwrite=True)\n",
" print(\"--->\", outfile)\n",
"# for i, (model, result) in enumerate(results.items()):\n",
"# result[1].data = flux_conv[:, None].value * (result[1].data - bgresult[1].data) / exptime\n",
"# result[1].header['BUNIT'] = \"e/s\"\n",
"# outfile = (Path(fitsfiles[model]).stem + \"-simulated_LSS_L-bgsub_fluxcal\" \n",
"# + Path(fitsfiles[model]).suffix)\n",
"# result.writeto(outfile, overwrite=True)\n",
"# print(\"--->\", outfile)\n",
"\n",
"for i, (model, result) in enumerate(rotresults.items()):\n",
" result[1].data = flux_conv[:, None].value * (result[1].data - bgresult[1].data) / exptime\n",
" result[1].header['BUNIT'] = \"e/s\"\n",
" outfile = (Path(fitsfiles[model]).stem + \"-rot_\" + str(angle) \n",
" + \"-simulated_LSS_L-bgsub_fluxcal\" + Path(fitsfiles[model]).suffix)\n",
" result.writeto(outfile, overwrite=True)\n",
" print(\"--->\", outfile)"
"# for i, (model, result) in enumerate(rotresults.items()):\n",
"# result[1].data = flux_conv[:, None].value * (result[1].data - bgresult[1].data) / exptime\n",
"# result[1].header['BUNIT'] = \"e/s\"\n",
"# outfile = (Path(fitsfiles[model]).stem + \"-rot_\" + str(angle) \n",
"# + \"-simulated_LSS_L-bgsub_fluxcal\" + Path(fitsfiles[model]).suffix)\n",
"# result.writeto(outfile, overwrite=True)\n",
"# print(\"--->\", outfile)"
]
},
{
Expand Down Expand Up @@ -764,12 +764,19 @@
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(12, 6))\n",
"plt.plot(det_wave, result_extract, label=\"flux-calibrated simulated spectrum\")\n",
"plt.plot(src_wave, src_extract, label=\"input spectrum\")\n",
"plt.legend()\n",
"plt.xlim(3.1, 4.0);"
"# plt.figure(figsize=(12, 6))\n",
"# plt.plot(det_wave, result_extract, label=\"flux-calibrated simulated spectrum\")\n",
"# plt.plot(src_wave, src_extract, label=\"input spectrum\")\n",
"# plt.legend()\n",
"# plt.xlim(3.1, 4.0);"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
Expand Down

0 comments on commit 4dd3b1b

Please sign in to comment.