#!/usr/bin/env python from argparse import ArgumentParser from argparse import RawTextHelpFormatter import numpy as np import math from scipy.ndimage.interpolation import rotate from astropy.io import fits import os parser=ArgumentParser( description='A script to produce datacubes and collapsed images from MocassinPlot outputs.\nInputs come from a tab-sepatated file named mocplot.in which must be placed in the input directory. An example is reproduced below\n\nsymmetricXYZ\tyes\npixels\t500\ndistance\t2.4\nincl\t15\n1\tHalpha\n2\tHbeta\n3\tHeI6678\n4\tHeII4686\n5\tCII4267\n6\tNII6583\n7\tOII3727\n8\tOII7332\n9\tOII7331\n10\tOIII5007\n11\tOIII4959\n12\tOIII4363\n\n\nsymmetricXYZ = yes or no. This defines whether the output is a single nebula or whether the mocassinPlot output should be considered a single quadrant and then grown into a full datacube. If not supplied, the model will be reproduced "as is".\n\npixels = int. The number of pixels per side in the newly formed data cube (N.B. the run time scales approximately exponentially with pixels!).\n\ndistance = float. Distance to the nebula in kpc, in order to place image on an arcsec scale. If not supplied, the default canonical value of 1 kpc will be used.\n\nincl = float. Inclination of the produced nebula in degrees c.f. mocassin grid. If not supplied, the model will not be inclined. Note, changing inclination will change the size and/or shape of the final datacube in order to contain all data.\n\nThe subsequent entries correspond to the outputs from mocassinPlot. In this example, plot.out contains data for twelve lines, the first (1) being in Hydrogen Alpha (assigned the label Halpha) and the last (12) OIII 4634A (assigned the label OIII4634). The assigned labels (the second column) are used by this program to name the resulting outputs. In order for the assigned labels to correspond to the correct outputs, the numeric index (the first column) must correspond to the order the lines were supplied to mocassinPlot.', epilog="A panther may love the ape, yet despite their efforts they will never produce issue.\nDavid Jones, 20 February 2015. djones@iac.es", formatter_class=RawTextHelpFormatter) args=parser.parse_args() kpc=3.08567758e21 #kpc in cm #Import parameters from mocplot.in try: with open("input/mocplot.in") as paramFile: params = [line.split() for line in paramFile] except: print "Couldn't find mocplot.in in the input folder!" quit() numplots=0 angle=0 ploti=1 dist=1 others=0 sym='no' for x in range (0,len(params)): if params[x][0] == 'symmetricXYZ': others=others+1 elif params[x][0] == 'distance': others=others+1 elif params[x][0] == 'pixels': others=others+1 elif params[x][0] == 'incl': others=others+1 numplots=len(params)-others plotnames=['blooper']*numplots for x in range (0,len(params)): if params[x][0] == 'symmetricXYZ': sym=params[x][1] elif params[x][0] == 'distance': dist=float(params[x][1]) elif params[x][0] == 'pixels': pix=int(params[x][1]) size=int(params[x][1]) elif params[x][0] == 'incl': angle=float(params[x][1]) else: try: params[x][0]=int(params[x][0]) except: print "Unrecognised input \"%s \t %s\", carrying on regardless..." % params[x][0],params[x][1] if params[x][0]==ploti: plotnames[ploti-1]=params[x][1] ploti=ploti+1 if size: data_cube = np.zeros ((size,size,size)) else: print "Size of required cube is undefined in mocplot.in" quit() print "Scaling axes for a distance of %r kpc" %dist try: with open("output/plot.out") as plotFile: plotcell = [[float(entry) for entry in line.split()] for line in plotFile] with open("output/grid4.out") as gridFile: gridcell = [[float(entry) for entry in line.split()] for line in gridFile] except: print "I couldn't find the outputs of mocassinPlot. Are you sure you've run it?" quit() #Find scale by looking for largest grid extent. Fits pixels will all be initialised to zero so grids that aren't a cube will just sit in the corner of the final fits cube. gridsize = len(gridcell) maxsize=0 for x in range (0,3): if gridcell[gridsize-1][x] > maxsize: maxsize = float(gridcell[gridsize-1][x])+gridcell[gridsize-1][3]**(1/3.0) maxarcsize = math.degrees(math.atan(maxsize/(dist*kpc)))*3600 if sym=='yes': print "Size of the output grid is %r cm or %r arcsec" %(2*maxsize,round(2*maxarcsize,3)) print "Mapping to all 8 quadrants" size=int(size/2) data_cube=np.zeros((numplots,2*size,2*size,2*size)) totalcells=(2*size)**3 plotcellvolume=(maxsize**3)/totalcells plotcellsizecm=maxsize/size plotcellsizearc = maxarcsize/size print "Total number of pixels to populate %s" %totalcells print "I will inform you periodically of my progress. If I go a long time without saying anything or returning the command line, something has probably gone awry." cellnum=0 rownumdummy=0 percent=5 for i in range (0,size): for j in range (0,size): for k in range (0,size): # Select best grid cell for this plot cell, attribute brightness mulitplied by plotcellvolume and divided by gridcellvolume brightness=0 cmx=i*plotcellsizecm cmy=j*plotcellsizecm cmz=k*plotcellsizecm for x in range (rownumdummy,gridsize): diffx=cmx-gridcell[x][0] rdiffx=round(diffx*10**(-15),3) rgcsize = round((gridcell[x][3])**(1/3.0),3) if abs(rdiffx) <= 1.1*rgcsize: diffy=cmy-gridcell[x][1] rdiffy=round(diffy*10**(-15),3) if abs(rdiffy) <= 1.1*rgcsize: diffz=cmz-gridcell[x][2] rdiffz=round(diffz*10**(-15),3) if abs(rdiffz) <= 1.1*rgcsize: cellsize=gridcell[x][3]*10**45 cellnum=cellnum+1 if x >= int(round(gridsize**(1/3.0),0)): rownumdummy=x-(size/2)*int(round(gridsize**(1/3.0),0)) if (100*8*cellnum/float(totalcells))>=percent: print "%i%% of grid filled!" %(int(100*8*cellnum/float(totalcells))) percent=percent+5 for y in range (0,numplots): data_cube[y][size+j][size+k][size+i]=plotcell[x][y+1]*plotcellvolume/cellsize data_cube[y][size-j][size+k][size+i]=plotcell[x][y+1]*plotcellvolume/cellsize data_cube[y][size+j][size-k][size+i]=plotcell[x][y+1]*plotcellvolume/cellsize data_cube[y][size+j][size+k][size-i]=plotcell[x][y+1]*plotcellvolume/cellsize data_cube[y][size-j][size-k][size+i]=plotcell[x][y+1]*plotcellvolume/cellsize data_cube[y][size+j][size-k][size-i]=plotcell[x][y+1]*plotcellvolume/cellsize data_cube[y][size-j][size+k][size-i]=plotcell[x][y+1]*plotcellvolume/cellsize data_cube[y][size-j][size-k][size-i]=plotcell[x][y+1]*plotcellvolume/cellsize break #Add "fast mode" where cells are all same size and map directly to number of pixels in final image. gridline=0, change for to while and add a gridline++1 at end of while loop else: print "Size of the output grid is %r cm or %r arcsec" %(maxsize,round(maxarcsize,3)) print "Not mapping to multiple quadrants" data_cube=np.zeros((numplots,size,size,size)) plotcellvolume=(maxsize**3)/(size**3) plotcellsizecm=maxsize/size plotcellsizearc = maxarcsize/size totalcells=size**3 print "Total number of pixels to populate %s" %totalcells print "I will inform you periodically of my progress. If I go a long time without saying anything or returning the command line, something has probably gone awry." cellnum=0 rownumdummy=0 percent=5 for i in range (0,size): for j in range (0,size): for k in range (0,size): # Select best grid cell for this plot cell, attribute brightness mulitplied by plotcellvolume and divided by gridcellvolume brightness=0 cmx=i*plotcellsizecm cmy=j*plotcellsizecm cmz=k*plotcellsizecm for x in range (rownumdummy,gridsize): diffx=cmx-gridcell[x][0] rdiffx=round(diffx*10**(-15),3) rgcsize = round((gridcell[x][3])**(1/3.0),3) if abs(rdiffx) <= 1.1*rgcsize: diffy=cmy-gridcell[x][1] rdiffy=round(diffy*10**(-15),3) if abs(rdiffy) <= 1.1*rgcsize: diffz=cmz-gridcell[x][2] rdiffz=round(diffz*10**(-15),3) if abs(rdiffz) <= 1.1*rgcsize: cellsize=gridcell[x][3]*10**45 cellnum=cellnum+1 if x >= int(round(gridsize**(1/3.0),0)): rownumdummy=x-(size/2)*int(round(gridsize**(1/3.0),0)) if (100*cellnum/float(totalcells))>=percent: print "%i%% of grid filled!" %(int(100*cellnum/float(totalcells))) percent=percent+5 for y in range (0,numplots): data_cube[y][j][k][i]=plotcell[x][y+1]*plotcellvolume/cellsize break #Print to datacubes try: os.system("rm moc_*.fits") os.system("rm output/moc_*.fits") except: pass for z in range (0,numplots): linedata_cube=data_cube[z] if angle: unrot_cube=linedata_cube linedata_cube=rotate(unrot_cube,angle,axes=(1,0),reshape=True) hdu = fits.PrimaryHDU(linedata_cube) hdu.header['CRPIX1']='0' hdu.header['CRPIX2']='0' hdu.header['CRPIX3']='0' hdu.header['CRVAL1']='0' hdu.header['CRVAL2']='0' hdu.header['CRVAL3']='0' hdu.header['CD1_1']='%r'%plotcellsizearc hdu.header['CD1_2']='0' hdu.header['CD1_3']='0' hdu.header['CD2_1']='0' hdu.header['CD2_2']='%r'%plotcellsizearc hdu.header['CD2_3']='0' hdu.header['CD3_1']='0' hdu.header['CD3_2']='0' hdu.header['CD3_3']='%r'%plotcellsizearc hdu.header['CTYPE1']='Angular size' hdu.header['CTYPE2']='Angular size' hdu.header['CTYPE3']='Angular size' hdu.writeto('moc_%scube.fits'%plotnames[z]) col_cube=linedata_cube[0] for x in range (0,pix): col_cube+=linedata_cube[x] hducol = fits.PrimaryHDU(col_cube) hducol.header['CRPIX1']='0' hducol.header['CRPIX2']='0' hducol.header['CD1_1']='%r'%plotcellsizearc hducol.header['CD1_2']='0' hducol.header['CD2_1']='0' hducol.header['CD2_2']='%r'%plotcellsizearc hducol.header['CRVAL1']='0' hducol.header['CRVAL2']='0' hducol.header['CTYPE1']='Angular size' hducol.header['CTYPE2']='Angular size' hducol.writeto('moc_%scol.fits'%plotnames[z]) os.system("mv moc_*.fits output/.")