library(astro)
## Loading required package: MASS
## Loading required package: plotrix
library(celestial)
## Loading required package: RANN
## Loading required package: NISTunits
## Loading required package: pracma
## 
## Attaching package: 'pracma'
## The following objects are masked from 'package:astro':
## 
##     eta, kron
library(ProFit)
## Loading required package: Rcpp
## Loading required package: magicaxis
## Loading required package: sm
## Package 'sm', version 2.2-5.4: type help(sm) for summary information
## 
## Attaching package: 'sm'
## The following object is masked from 'package:pracma':
## 
##     nile
## The following object is masked from 'package:MASS':
## 
##     muscle
## Loading required package: mapproj
## Loading required package: maps
## Loading required package: FITSio
## Loading required package: data.table
matt_image=read.fits('~/Downloads/9008500074_vst_400.fits')

Have a look (we will focus on the r-band data for this example). Note we need to subtract 0.5 because R uses half integer pixel centres and FITS uses integer pixel centres. This means when converting back to RA/DEC from ProFit xy we will need to add 0.5.

magimageWCS(matt_image$dat[[5]], matt_image$hdr[[5]])
sami_ra=as.numeric(get.fitskey('SAMIRA',matt_image$hdr[[5]]))
sami_dec=as.numeric(get.fitskey('SAMIDEC',matt_image$hdr[[5]]))
sami_xy=radec2xy(sami_ra,sami_dec,header = matt_image$hdr[[5]])-0.5
points(sami_xy[1],sami_xy[2],col='red', pch=3)

The object of interest is pretty much exactly in the centre of the image.

Also worth checking the weight map:

magimageWCS(matt_image$dat[[6]], matt_image$hdr[[5]])

There is a large variation in the stack size throughout the frame, especially in the region of our galaxy of interest. This means it is very important we use the weight map (actually an inverse sky variacen map) when consructing our sigma map.

Do our basic segmentation using the inverse variance weight map provided:

matt_profound=profitProFound(matt_image$dat[[5]], tolerance = 10, sigma=3, magzero=0, pixscale=0.2, sky=0, skyRMS=sqrt(1/matt_image$dat[[6]]), plot=TRUE, verbose=TRUE, redosky=FALSE)
## [1] "Running profitProFound:"
## [1] "Skipping making initial sky map - User provided sky and sky RMS"
## [1] "Making initial segmentation image - 0 sec"
## [1] " - Running profitMakeSegim:"
## [1] " - Skipping making initial local estimate of the sky - User provided sky"
## [1] " - Skipping making initial local estimate of the sky RMS - User provided sky RMS"
## [1] " - Smoothing the image - 0.409999999999999 sec"
## [1] " - Watershed de-blending - 0.742 sec"
## [1] " - Skipping segmentation plot - plot set to FALSE"
## [1] " - Skipping making final local estimate of the sky - User provided sky"
## [1] " - Skipping making initial local estimate of the sky RMS - User provided sky RMS"
## [1] " - Skipping segmentation statistics - segstats set to FALSE"
## [1] " - profitMakeSegim is finished! - 1.296 sec"
## [1] "Skipping making better sky map - User provided sky and sky RMS"
## [1] "Calculating initial segstats - 1.296 sec"
## [1] "Doing dilations:"
## [1] "Iteration 1 of 6 - 1.798 sec"
## [1] " - Running profitMakeSegimDilate:"
## [1] " - Dilating segments - 0.000999999999999446 sec"
## [1] " - Calculating segstats - 0.841 sec"
## [1] " - profitMakeSegimDilate is finished! - 1.285 sec"
## [1] "Iteration 2 of 6 - 3.11 sec"
## [1] " - Running profitMakeSegimDilate:"
## [1] " - Dilating segments - 0.000999999999999446 sec"
## [1] " - Calculating segstats - 0.753 sec"
## [1] " - profitMakeSegimDilate is finished! - 1.167 sec"
## [1] "Iteration 3 of 6 - 4.331 sec"
## [1] " - Running profitMakeSegimDilate:"
## [1] " - Dilating segments - 0 sec"
## [1] " - Calculating segstats - 0.757999999999999 sec"
## [1] " - profitMakeSegimDilate is finished! - 1.135 sec"
## [1] "Iteration 4 of 6 - 5.493 sec"
## [1] " - Running profitMakeSegimDilate:"
## [1] " - Dilating segments - 0 sec"
## [1] " - Calculating segstats - 0.930000000000001 sec"
## [1] " - profitMakeSegimDilate is finished! - 1.315 sec"
## [1] "Iteration 5 of 6 - 6.837 sec"
## [1] " - Running profitMakeSegimDilate:"
## [1] " - Dilating segments - 0 sec"
## [1] " - Calculating segstats - 0.811 sec"
## [1] " - profitMakeSegimDilate is finished! - 1.167 sec"
## [1] "Iteration 6 of 6 - 8.032 sec"
## [1] " - Running profitMakeSegimDilate:"
## [1] " - Dilating segments - 0 sec"
## [1] " - Calculating segstats - 0.850000000000001 sec"
## [1] " - profitMakeSegimDilate is finished! - 1.261 sec"
## [1] "Finding CoG convergence - 9.321 sec"
## [1] "Constructing final segim - 9.34 sec"
## [1] "Skipping making final sky map - redosky set to FALSE"
## [1] "Calculating final segstats - 10.421 sec"
## [1] "Plotting segments - 10.748 sec"

## [1] "profitProFound is finished! - 18.312 sec"

Looks to have worked, making proper use of the sky weights. We can check the distribution of surface brightness limits:

maghist(matt_profound$SBlim)
## [1] "Summary of used sample:"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   22.52   23.01   23.01   22.93   23.02   23.02 
## [1] "sd / 1-sig / 2-sig range:"
## [1] 0.1529431 0.1878388 0.1929730
## [1] "Plotting 3972293 out of 4004001 (99.21%) data points (0 < xlo & 0 > xhi)"

Mostly at a 2-skyRMS depth of 23 mag/asec^2.

Next we want to grab the gain from the header:

matt_gain=as.numeric(matt_image$hdr[[5]][matt_image$hdr[[5]][,'key']=='GAIN','value'])
matt_gain
## [1] 1.5248e+12

Now we make our sigma map:

matt_sigma=profitMakeSigma(matt_image$dat[[5]], object=matt_profound$objects, sky=0, skyRMS=sqrt(1/matt_image$dat[[6]]), gain=matt_gain, plot = TRUE)

That is one funky sigma map! Anyway, it all looks pretty legit.

Now let’s select some stars:

magplot(matt_profound$segstats$R50, matt_profound$segstats$SB_N90, log='x', col=hsv(magmap(matt_profound$segstats$axrat, flip=TRUE)$map), xlab='Major Axis R50 / asec', ylab='Surface Brightness within R90 mag/asec^2')
magbar('topright', title='Axrat', titleshift=1)
abline(h=22, lty=2)
abline(v=c(0.75,0.9), lty=2)

Wow, nice data. I spy some stars. That’s probably a good enough selection for now.

matt_starlist=matt_profound$segstats[matt_profound$segstats$R50>0.75 & matt_profound$segstats$R50<0.9 & matt_profound$segstats$SB_N90<22 & matt_profound$segstats$axrat>0.9,]
matt_starlist
##    segID      xcen      ycen         flux      mag    N N50 N90        R
## 1      1 1333.4249  558.5962 6.397906e-06 12.98491 5879  59 556 8.822063
## 2      2 1309.3200  788.0734 2.720483e-06 13.91339 3744  58 497 7.038853
## 3      3  422.7748 1034.9577 1.091321e-06 14.90512 2269  56 438 5.453776
## 4      4 1722.7194 1799.5817 4.976788e-07 15.75763 1315  52 353 4.221532
## 5      5 1362.6744 1236.5223 3.476700e-07 16.14708 1317  55 365 4.179721
## 7      7  952.6647  155.2395 1.902370e-07 16.80176 1053  55 346 3.747293
## 10    10  990.1668  731.6007 2.065516e-07 16.71243  956  56 316 3.516234
## 11    11 1459.0088  519.4833 1.357174e-07 17.16841  909  52 313 3.458318
## 13    13 1575.6473  404.9535 1.252445e-07 17.25560  884  53 322 3.425370
## 19    19 1607.8759  420.2830 6.784138e-08 17.92126  696  52 291 3.033739
## 24    24  504.5669 1563.6473 3.990464e-08 18.49744  543  49 242 2.695103
## 25    25 1698.2141 1733.5040 3.674489e-08 18.58701  514  48 221 2.639632
## 26    26  736.0374 1689.8408 3.589143e-08 18.61252  516  50 236 2.619551
## 28    28  177.9207 1537.4079 3.245027e-08 18.72195  496  49 232 2.562273
## 33    33 1129.9299 1272.7250 2.314022e-08 19.08908  656  44 196 2.952206
## 35    35  172.2343 1770.8807 1.625648e-08 19.47243  481  47 202 2.497374
## 42    42 1309.5631  393.8281 1.146962e-08 19.85113  407  44 162 2.343286
## 52    52 1114.6144 1105.8782 1.240191e-08 19.76628  499  43 162 2.601340
##          R50      R90     SB_N   SB_N50   SB_N90      xsd      ysd
## 1  0.8837805 2.713038 18.91331 14.66976 16.46714 2.167403 2.188411
## 2  0.8760885 2.564556 19.35187 15.57968 17.27382 2.178220 2.203580
## 3  0.8567893 2.396167 19.79985 16.53331 18.12835 2.156133 2.182541
## 4  0.8394772 2.187231 20.06009 17.30536 18.74661 2.140435 2.119404
## 5  0.8541532 2.200396 20.45120 17.75571 19.17236 2.191926 2.220013
## 7  0.8564161 2.148036 20.86298 18.41039 19.76900 2.267334 2.233567
## 10 0.8510259 2.021586 20.66872 18.34062 19.58119 2.269347 2.300306
## 11 0.8271507 2.029342 21.06997 18.71614 20.02682 2.198637 2.197648
## 13 0.8387245 2.067328 21.12688 18.82402 20.14479 2.197767 2.212996
## 19 0.8292304 1.961643 21.53294 19.46900 20.70054 2.245549 2.236395
## 24 0.8096057 1.799216 21.83959 19.98066 21.07652 2.177556 2.265241
## 25 0.8066451 1.730844 21.86957 20.04784 21.06753 2.217459 2.249391
## 26 0.8154307 1.771569 21.89930 20.11767 21.16435 2.243959 2.293256
## 28 0.8053461 1.752381 21.96581 20.20517 21.25522 2.230894 2.239853
## 33 0.7645768 1.613700 22.63649 20.45544 21.43927 2.471520 2.433613
## 35 0.7806566 1.618403 22.68295 20.91040 21.85536 2.394471 2.419594
## 42 0.7704678 1.478379 22.88026 21.21748 21.99446 2.495590 2.355425
## 52 0.7636267 1.482192 23.01668 21.10767 21.90961 2.598527 2.651808
##          covxy        corxy       con asymm flux_reflect mag_reflect
## 1   0.17902718  0.037744204 0.3257531    NA           NA          NA
## 2   0.17659688  0.036791913 0.3416141    NA           NA          NA
## 3   0.12449287  0.026454920 0.3575666    NA           NA          NA
## 4   0.27919604  0.061545095 0.3838083    NA           NA          NA
## 5   0.18945920  0.038934472 0.3881816    NA           NA          NA
## 7  -0.22153442 -0.043744822 0.3986973    NA           NA          NA
## 10  0.04064985  0.007787040 0.4209693    NA           NA          NA
## 11  0.15853467  0.032810484 0.4075956    NA           NA          NA
## 13  0.19924902  0.040966972 0.4057047    NA           NA          NA
## 19  0.18886754  0.037608510 0.4227223    NA           NA          NA
## 24  0.14614877  0.029628625 0.4499770    NA           NA          NA
## 25  0.30392015  0.060931089 0.4660414    NA           NA          NA
## 26  0.19383854  0.037668000 0.4602873    NA           NA          NA
## 28  0.19284068  0.038592259 0.4595725    NA           NA          NA
## 33  0.23831359  0.039621718 0.4738035    NA           NA          NA
## 35 -0.08650302 -0.014930667 0.4823622    NA           NA          NA
## 42  0.01972289  0.003355274 0.5211573    NA           NA          NA
## 52  0.41084349  0.059622080 0.5152010    NA           NA          NA
##         maj      min     axrat       ang iter  origfrac
## 1  2.219949 2.135090 0.9617745 142.16794    2 0.9980113
## 2  2.232786 2.148272 0.9621486 143.73221    2 0.9964849
## 3  2.200736 2.137559 0.9712926 147.35502    2 0.9939067
## 4  2.195321 2.062499 0.9394976 130.44258    2 0.9894373
## 5  2.250740 2.160362 0.9598454 144.05443    2 0.9897662
## 7  2.301960 2.197864 0.9547797  54.46643    2 0.9849511
## 10 2.302663 2.266955 0.9844929 165.05787    2 0.9846776
## 11 2.233916 2.161777 0.9677076 134.60719    2 0.9811486
## 13 2.250739 2.159098 0.9592842 139.78404    2 0.9794646
## 19 2.282971 2.198181 0.9628598 131.90066    2 0.9678486
## 24 2.275972 2.166337 0.9518293 161.55933    2 0.9601703
## 25 2.302307 2.162468 0.9392612 141.60303    2 0.9633485
## 26 2.317537 2.218874 0.9574278 149.99151    2 0.9598637
## 28 2.278331 2.191583 0.9619248 137.96466    2 0.9538734
## 33 2.504246 2.399924 0.9583420 124.34428    4 0.8713944
## 35 2.428890 2.385041 0.9819467  27.52193    3 0.8613470
## 42 2.495704 2.355304 0.9437432  91.66011    3 0.8573074
## 52 2.706698 2.541302 0.9388936 144.40053    4 0.7698467

Take a look at some stars:

for (i in 1:length(matt_starlist[,1])){
  magimage(magcutout(matt_image$dat[[5]], loc=matt_starlist[i,c('xcen','ycen')])$image)
}

In practice you might want to fit all of these, but we will take the fifth (not saturated):

psf_image=magcutout(matt_image$dat[[5]], loc=matt_starlist[5,c('xcen','ycen')])
psf_sigma=magcutout(matt_sigma, loc=matt_starlist[5,c('xcen','ycen')])
psf_segim=magcutout(matt_profound$segim, loc=matt_starlist[5,c('xcen','ycen')])

Look at our target star to check we are happy:

magimage(psf_image$image)

magimage(psf_sigma$image)

magimage(psf_segim$image)

Next we make our initial model, where we take our PSF estimates from the segmentation stats:

psf_x=psf_image$loc[1]
psf_y=psf_image$loc[2]
psf_mag=matt_starlist[5,'mag']
psf_fwhm=matt_starlist[5,'R50']*2/0.2
psf_con=1/matt_starlist[5,'con']

psf_modellist=list(
  moffat=list(
    xcen=psf_x,
    ycen=psf_y,
    mag=psf_mag,
    fwhm=psf_fwhm,
    con=psf_con,
    axrat=1,
    box=0
  )
)

We also need to set up the tofit, tolog and intervals structure

psf_tofit=list(
  moffat=list(
    xcen=TRUE,
    ycen=TRUE,
    mag=TRUE,
    fwhm=TRUE,
    con=TRUE,
    axrat=FALSE,
    box=FALSE
  )
)

psf_tolog=list(
  moffat=list(
    xcen=FALSE,
    ycen=FALSE,
    mag=FALSE,
    fwhm=TRUE,
    con=TRUE,
    axrat=FALSE,
    box=FALSE
  )
)

psf_intervals=list(
  moffat=list(
    xcen=list(psf_x+c(-5,5)),
    ycen=list(psf_y+c(-5,5)),
    mag=list(c(psf_mag-2,psf_mag+2)),
    fwhm=list(c(1,100)),
    con=list(c(1,100)),
    axrat=list(c(0.1,1)),
    box=list(c(-1,1))
  )
)

Setup our fitting structure:

psf_Data=profitSetupData(psf_image$image, sigma=psf_sigma$image, modellist=psf_modellist, tofit=psf_tofit, tolog=psf_tolog, magzero=0, algo.func='optim', intervals=psf_intervals, segim=psf_segim$image)

We can check our intial guess easily:

profitLikeModel(parm=psf_Data$init, Data=psf_Data, makeplots=TRUE, plotchisq=TRUE)
## Warning in xy.coords(x, y, xlabel, ylabel, log): 13 y values <= 0 omitted
## from logarithmic plot

## [1] -3993.263

Not bad, but not perfect. We can do an optim fit to improve things:

psf_fit=optim(psf_Data$init, profitLikeModel, method='BFGS', Data=psf_Data, control=list(fnscale=-1))

Now we plot the output:

profitLikeModel(parm=psf_fit$par, Data=psf_Data, makeplots=TRUE, plotchisq=TRUE)
## Warning in xy.coords(x, y, xlabel, ylabel, log): 14 y values <= 0 omitted
## from logarithmic plot
## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 y value <= 0 omitted
## from logarithmic plot

## [1] -1941.373

This looks like a very good fit (hint of ellipticity, so might want to fit for that too). We will now contruct our model PSF:

psf_modellist_fit=profitRemakeModellist(parm=psf_fit$par, Data=psf_Data)$modellist
psf_modellist_fit$moffat$xcen=25/2
psf_modellist_fit$moffat$ycen=25/2
psf_model=profitMakeModel(modellist=psf_modellist_fit, magzero=30, dim=c(25,25))$z

And now we can look at our final PSF model:

magimage(psf_model)

Note the FWHM is 5.7 pixels, which is 1.14 asec, which sounds about right for VST data.

Now we can cutout our galaxy of interest, and expand the segmentation map slightly:

galrow=which.min(sqrt((matt_profound$segstats$xcen-dim(matt_image$dat[[5]])[1]/2)^2+(matt_profound$segstats$ycen-dim(matt_image$dat[[5]])[2]/2)^2))
matt_profound$segstats[galrow,]
##    segID     xcen     ycen         flux      mag     N  N50  N90        R
## 14    14 997.8289 1000.319 1.336834e-06 14.68481 10880 1797 6587 15.33741
##         R50      R90     SB_N  SB_N50   SB_N90      xsd     ysd     covxy
## 14 6.233209 11.93387 21.28153 20.0789 20.85107 9.792726 16.5685 -9.633415
##          corxy       con asymm flux_reflect mag_reflect      maj      min
## 14 -0.05937363 0.5223122    NA           NA          NA 16.58412 9.766239
##        axrat      ang iter origfrac
## 14 0.5888909 3.078242    2 0.981575

Much like we did with the PSF fit, we have to extract the region of interest a setup our fitting structures:

gal_image=magcutout(matt_image$dat[[5]], loc=matt_profound$segstats[galrow,c('xcen','ycen')], box=c(301,301))
gal_sigma=magcutout(matt_sigma, loc=matt_profound$segstats[galrow,c('xcen','ycen')], box=c(301,301))
gal_segim=magcutout(matt_profound$segim, loc=matt_profound$segstats[galrow,c('xcen','ycen')], box=c(301,301))

Check the cutout regions:

magimage(gal_image$image)

magimage(gal_sigma$image)

magimage(gal_segim$image)

The sky looks to be a bit high TBH. I will subtract the ProFit estimate:

gal_image$image=gal_image$image-profitSkyEst(gal_image$image, objects=gal_segim$image>0,plot=TRUE)$sky

Now we expand out the region for fitting:

gal_segim_expand=profitMakeSegimExpand(gal_image$image, segim=gal_segim$image, expand=matt_profound$segstats[galrow,'segID'], expandsigma=20, skycut=-2, magzero=0, pixscale=0.2, sky=0, plot=TRUE)

## [1] " - Calculating segstats - 0.623000000000005 sec"

We will select the object from the subset (usually this will be the first entry, but it might not be if there is a very bright sources nearby):

galrow_subset=which.min(sqrt((gal_segim_expand$segstats$xcen-gal_image$loc[1])^2+(gal_segim_expand$segstats$ycen-gal_image$loc[2])^2))
gal_segim_expand$segstats[galrow_subset,]
##   segID     xcen     ycen         flux      mag     N  N50   N90        R
## 1    14 149.9939 150.7263 1.564792e-06 14.51386 39312 2338 12299 27.32732
##        R50      R90     SB_N   SB_N50   SB_N90     xsd      ysd     covxy
## 1 6.664333 15.28514 22.50532 20.19369 21.35808 12.6012 18.74668 -11.00842
##         corxy       con asymm flux_reflect mag_reflect     maj      min
## 1 -0.04660029 0.4360009    NA           NA          NA 18.7634 12.57629
##       axrat      ang
## 1 0.6702568 3.259891

As before, we have to setup some reasonable guesses for our galaxy fit:

gal_x=gal_image$loc[1]
gal_y=gal_image$loc[2]
gal_mag=gal_segim_expand$segstats[galrow_subset,'mag']
gal_re=gal_segim_expand$segstats[galrow_subset,'R50']/0.2
gal_ang=gal_segim_expand$segstats[galrow_subset,'ang']
gal_axrat=gal_segim_expand$segstats[galrow_subset,'axrat']

gal_modellist=list(
  sersic=list(
    xcen=rep(gal_x,2),
    ycen=rep(gal_y,2),
    mag=c(gal_mag+4, gal_mag),
    re=c(gal_re/4,gal_re),
    nser=c(4,1),
    ang=c(0,gal_ang),
    axrat=c(1,gal_axrat),
    box=rep(0,2)
  )
)

And now we make the other inputs we need for our fitting Data structure:

gal_tofit=list(
  sersic=list(
    xcen= c(TRUE,NA), #We fit for xcen and tie the two togther
    ycen= c(TRUE,NA), #We fit for ycen and tie the two togther
    mag= c(TRUE,TRUE), #Fit for both
    re= c(TRUE,TRUE), #Fit for both
    nser= c(TRUE,FALSE), #Fit for bulge
    ang= c(FALSE,TRUE), #Fit for disk
    axrat= c(FALSE,TRUE), #Fit for disk
    box= c(FALSE,FALSE) #Fit for neither
  )
)

gal_tolog=list(
  sersic=list(
    xcen= c(FALSE,FALSE),
    ycen= c(FALSE,FALSE),
    mag= c(FALSE,FALSE),
    re= c(TRUE,TRUE), #re is best fit in log space
    nser= c(TRUE,TRUE), #nser is best fit in log space
    ang= c(FALSE,FALSE),
    axrat= c(TRUE,TRUE), #axrat is best fit in log space
    box= c(FALSE,FALSE)
  )
)

gal_intervals=list(
  sersic=list(
    xcen=list(lim=gal_x+c(-5,5),lim=gal_x+c(-5,5)),
    ycen=list(lim=gal_y+c(-5,5),gal_y+c(-5,5)),
    mag=list(lim=c(10,30),lim=c(10,30)),
    re=list(lim=c(0.1,gal_re),lim=c(0.1,gal_re*2)),
    nser=list(lim=c(0.5,20),lim=c(0.5,20)),
    ang=list(lim=c(-180,360),lim=c(-180,360)),
    axrat=list(lim=c(0.001,1),lim=c(0.001,1)),
    box=list(lim=c(-1,1),lim=c(-1,1))
  )
)

Now we can set up our full fitting structure with averything we have made:

gal_Data=profitSetupData(gal_image$image, sigma=gal_sigma$image, modellist=gal_modellist, tofit=gal_tofit, tolog=gal_tolog, magzero=0, algo.func='optim', intervals=gal_intervals, psf=psf_model, segim=gal_segim_expand$segim)
profitLikeModel(parm=gal_Data$init, Data=gal_Data, makeplots=TRUE, plotchisq=TRUE)

## [1] -67644.28

Not a bad start, but certainly not great. Let’s do an optim fit:

We can now recreate our optimised modellist:

gal_fit_modellist=profitRemakeModellist(gal_fit$par, Data=gal_Data)$modellist
print(gal_fit_modellist)
## $sersic
## $sersic$xcen
## [1] 152.7955 152.7955
## 
## $sersic$ycen
## [1] 150.3614 150.3614
## 
## $sersic$mag
## [1] 17.41599 14.61961
## 
## $sersic$re
## [1]  1.530662 39.390160
## 
## $sersic$nser
## [1] 0.5 1.0
## 
## $sersic$ang
## [1] 0.000000 2.614332
## 
## $sersic$axrat
## [1] 1.0000000 0.5923211
## 
## $sersic$box
## [1] 0 0

And we can then plot the output:

profitLikeModel(parm=gal_fit$par, Data=gal_Data, makeplots=TRUE, plotchisq=TRUE)

## [1] -62297.63
profitEllipsePlot(Data=gal_Data, modellist=gal_fit_modellist, pixscale=0.2, SBlim=median(matt_profound$SBlim))

Not a bad first autoamtic effort. It looks liek you have some sky subtraction issues though (note a lot of red sky pixels in the bottom-right plot, there should be a simialr number of blue and red pixels for well subtracted sky).