use PDL::Graphics::Simple for plotting, eliminate loop

statocles-site/blog/2024/12/25/jwst/index.markdown

@@ -50,7 +50,7 @@ This is in the FITS format which is a RAW-type high dynamic range format used fo
 OK let's make a plot:
 
     use PDL;
-    use PDL::Graphics::Simple;
+    use PDL::Graphics::Simple qw(imag line plot);
     $jwst = rfits 'goodsn-wide0-v3_prism-clear_1211_1268.spec.fits[2]'; # Image is in second FITS extension
     imag $jwst, -0.05, 0.1, {justify=>0};
     @dims = dims($jwst);
@@ -80,7 +80,6 @@ We need to only extract the central few rows. Code like
 `$jwst->slice(':,13:17')->t()->sumover`  would work, but we can find a more optimal solution and demonstrate 'broadcasting' at the same time! First look at this code:
 
     # Make a gaussian extraction profile
-
     $x = xvals($dims[1]);
     $fwhm = 2.1; # pix
     $gaussian = 0.04*exp(-0.5*(($x-15)/($fwhm/2.35))**2);
@@ -120,26 +119,29 @@ We can now look at our resulting spectrum:
 
     # Get wavelength array from fits table
     $temp = rfits 'goodsn-wide0-v3_prism-clear_1211_1268.spec.fits[1]';
-    $wavelengths = $temp->{'wave    '};
-
-    # Make a nice plot
-
-    env min($wavelengths), max($wavelengths), -0.1, 1.1, {XTitle=>'Wavelength / \gmm', YTitle=>'Flux'};
-    line $wavelengths, $spec1d;
-    hold;
+    $wavelengths = $temp->{wave}; # if no Astro::FITS::Header, need 'wave    '
 
     @line_waves = (1216, 3728, 3870, 3890, 3971, 4103, 4342);  # In Angstroms
-    @line_names = ('Ly\ga','[OII]','[NeIII]', 'H\gz', 'H\ge', 'H\gd', 'H\gg');
+    # UTF-8, in bytes which is what Gnuplot needs
+    @line_names = ('|Lyα','|[OII]','|[NeIII]', '|Hζ', '|Hε', '|Hδ', '|Hγ');
+
+    $w_obs = pdl(\@line_waves)/10000; # Angstroms -> Microns
     $redshift = 10.603;
-    for $w (@line_waves){
-       $name = shift(@line_names);
-       $w_obs = ($w/10000) * (1+$redshift); # Angstroms -> Microns and apply the redshift.
-       line pdl($w_obs,$w_obs), pdl(-1,1), {Color=>Red};
-       text $name,$w_obs, 1.02, {Color=>Red, JUSTIFICATION=>0.5};
-    }
+    $w_obs *= 1+$redshift; # apply the redshift
+    $label_height = $w_obs->xlinvals(0.8,1.02);
+    my @heights = ($label_height - 0.02)->list;
+
+    # Make a nice plot
+    plot
+      with=>'lines', $wavelengths, $spec1d,
+      (map +(with=>'lines', style=>7, pdl($_,$_), pdl(-1,shift @heights)), $w_obs->list),
+      with=>'labels', style=>7, $w_obs, $label_height, \@line_names,
+      {xlabel=>'Wavelength / μm', ylabel=>'Flux',
+        xrange=>[minmax($wavelengths)],yrange=>[-0.1,1.1]};
+
+This code adds the wavelength calibration of the axis (which is complex and beyond the scope of this article) and makes a plot, adding labels for the location of some common spectral lines. Here is the plot:
 
-This code adds the wavelength calibration of the axis (which is complex and beyond the scope of this article) and makes a plot, adding labels for the location of some common spectral lines [OK I had to write a trivial loop for this :-(]. Here is the plot:
-![spectrum](spectrum.jpg)
+![spectrum](spectrum.png)
 
 What we see is a spectrum running from wavelengths 1 to 5 micrometers - this is beyond the range that the human eye can see because JWST is an infrared telescope. Looking at the spectral transitions there are two notable astrophysical features:
 
@@ -150,4 +152,4 @@ What we see is a spectrum running from wavelengths 1 to 5 micrometers - this is
 
 1. If you want to know more about this galaxy and its spectrum you can read the wikipedia page linked above or take a look at this (admittedly technical) [scientific paper](https://www.aanda.org/articles/aa/full_html/2023/09/aa46159-23/aa46159-23.html). You will see a similar spectrum in the paper.
 
-2. An interesting and straight forward exercise would be to use PDL to extract the slit profile and fit it with a gaussian to generation the extraction profile, rather than use my pre-canned one.
+2. An interesting and straight forward exercise would be to use PDL to extract the slit profile and fit it with a gaussian to generate the extraction profile, rather than use my pre-canned one.