diff --git a/.github/workflows/main.yaml b/.github/workflows/main.yaml index b025f11..cf9229c 100644 --- a/.github/workflows/main.yaml +++ b/.github/workflows/main.yaml @@ -6,10 +6,7 @@ on: jobs: build: - runs-on: ubuntu-latest - - steps: - uses: actions/checkout@v4 - name: Set up Python @@ -30,4 +27,6 @@ jobs: - name: Test with pytest run: | pip install pytest pytest-cov - pytest --doctest-modules --junitxml=junit/test-results.xml --cov=com --cov-report=xml --cov-report=html + pytest --junitxml=junit/test-results.xml + # pip install pytest-cov + # pytest --doctest-modules --junitxml=junit/test-results.xml --cov=com --cov-report=xml --cov-report=html diff --git a/.gitignore b/.gitignore index 6fcb2a5..82e1ae8 100644 --- a/.gitignore +++ b/.gitignore @@ -61,3 +61,4 @@ target/ .DS_Store **/.idea/* pyhdust/.* +junit/* diff --git a/docs/index.rst b/docs/index.rst index 41e98f0..2b6a396 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -4,7 +4,7 @@ Welcome to **pyhdust** documentation! ============================================================= **Analysis tools for multi-technique astronomical data and** *hdust* **models**. About the *hdust* code, see Carciofi & Bjorkman (`2006 `_, `2008 `_). -Pyhdust is currently at **version 1.5.12**. +Pyhdust is currently at **version 1.5.12-1**. **pyhdust** should be independent of plataform (Linux, Mac, Windows) and compatible with any version of Python (3.6+). diff --git a/scripts/extract3Dfits.py b/scripts/extract3Dfits.py index 8fe4f74..c4fade4 100755 --- a/scripts/extract3Dfits.py +++ b/scripts/extract3Dfits.py @@ -7,28 +7,30 @@ import sys from argparse import ArgumentParser import os -import pyfits as pf +import astropy.io.fits as pf -__version__ = "0.91" +__version__ = "0.92" __author__ = "Daniel Moser" __email__ = "dmfaes@gmail.com" -class MyParser(ArgumentParser): +class MyParser(ArgumentParser): def error(self, message): - sys.stderr.write('# ERROR! %s\n' % message) + sys.stderr.write("# ERROR! %s\n" % message) self.print_help() sys.exit(2) + parser = MyParser(description=__doc__) -parser.add_argument('--version', action='version', - version='%(prog)s {0}'.format(__version__)) +parser.add_argument( + "--version", action="version", version="%(prog)s {0}".format(__version__) +) parser.add_argument("INPUT", help=("Filename of the 3D FITS cube"), type=str) args = parser.parse_args() -if __name__ == '__main__': +if __name__ == "__main__": fname = args.INPUT fits = pf.open(fname) @@ -39,6 +41,7 @@ def error(self, message): hdulist[0].header = fits[0].header path, oname = os.path.split(fname) pref, ext = os.path.splitext(oname) - hdu.writeto(os.path.join(path, pref+"_{0:04d}".format(i)+ext), - overwrite=True) - print('# Saved {0} 2D images from {1}'.format(nimgs, oname)) + hdu.writeto( + os.path.join(path, pref + "_{0:04d}".format(i) + ext), overwrite=True + ) + print("# Saved {0} 2D images from {1}".format(nimgs, oname)) diff --git a/scripts/extract_ms_opera.py b/scripts/extract_ms_opera.py index 2ed8d9a..12a204c 100755 --- a/scripts/extract_ms_opera.py +++ b/scripts/extract_ms_opera.py @@ -13,12 +13,12 @@ import matplotlib.pyplot as plt from scipy.signal import savgol_filter -__version__ = "0.9" +__version__ = "0.91" __author__ = "Daniel Moser" __email__ = "dmfaes@gmail.com" -def plot_orders(ax, fwl, fflx, mode='o'): +def plot_orders(ax, fwl, fflx, mode="o"): for i in range(len(fwl)): ax.plot(fwl[i], fflx[i], mode, ms=1) return ax @@ -28,12 +28,12 @@ def do_norm(wl, flx): xout, yout = spt.max_func_pts(wl, flx, ws=0.01, avgbox=5) tck = interpolate.splrep(xout, yout) ynew = interpolate.splev(xout, tck, der=1) - dmed = np.median( np.abs(ynew) ) + dmed = np.median(np.abs(ynew)) l = len(wl) - # ConditionA: small derivative - conda = ( np.abs(ynew) < dmed ) + # ConditionA: small derivative + conda = np.abs(ynew) < dmed # ConditionB: in the limits - condb = ( wl[int(l*0.15)] > xout ) | ( xout > wl[int(l*(1-.15))] ) + condb = (wl[int(l * 0.15)] > xout) | (xout > wl[int(l * (1 - 0.15))]) # xmax = phc.find_nearest(yout, np.max(yout), idx=True) # xmax = xout[xmax] # xmax = wl[l/2] @@ -41,9 +41,9 @@ def do_norm(wl, flx): # condc = ( xout <= xmax ) & ( ynew >= 0 ) # ConditionD: xout > max(flx) AND ynew < 0 # condd = ( xout >= xmax ) & ( ynew <= 0 ) - conde = ( ( xout > 658.45 ) | ( xout < 654.1 ) ) - # idx = np.where( (( conda & (condd | condc) ) | condb) & (conde) ) - idx = np.where( (conda | condb) & (conde) ) + conde = (xout > 658.45) | (xout < 654.1) + # idx = np.where( (( conda & (condd | condc) ) | condb) & (conde) ) + idx = np.where((conda | condb) & (conde)) # & (condd | condc)) xout = xout[idx] yout = yout[idx] @@ -54,20 +54,20 @@ def do_norm(wl, flx): return wl, flx, ynew -path = '.' -opsps = glob('*.spc.gz') -print opsps +path = "." +opsps = glob("*.spc.gz") +print(opsps) avgspecs = [] for sp in opsps: - os.system('gunzip -c {0} > tmp.txt'.format(sp)) - spec = np.loadtxt('tmp.txt', skiprows=11) + os.system("gunzip -c {0} > tmp.txt".format(sp)) + spec = np.loadtxt("tmp.txt", skiprows=11) fig, ax = plt.subplots() - ax.plot(spec[:, 5], spec[:, 10], label='Leg.') + ax.plot(spec[:, 5], spec[:, 10], label="Leg.") # ax.set_title('Title') # ax.legend() - plt.savefig(sp.replace('.spc.gz', '_ori')) + plt.savefig(sp.replace(".spc.gz", "_ori")) plt.close(fig) # i = 0 @@ -83,46 +83,50 @@ def do_norm(wl, flx): try: wl, flx, ynew = do_norm(wl, flx) fwl.append(wl) - onorm.append(flx/ynew) + onorm.append(flx / ynew) ocont.append(ynew) - oflx.append(flx) + oflx.append(flx) except: - print('# Error in order {0}'.format(o)) + print("# Error in order {0}".format(o)) fig, ax = plt.subplots() - ax = plot_orders(ax, fwl, oflx, 'o') - ax = plot_orders(ax, fwl, ocont, '-') - plt.savefig(sp.replace('.spc.gz', '.orders')) + ax = plot_orders(ax, fwl, oflx, "o") + ax = plot_orders(ax, fwl, ocont, "-") + plt.savefig(sp.replace(".spc.gz", ".orders")) plt.close(fig) swl, sflx = spt.sum_ec(fwl, oflx) swl, scont = spt.sum_ec(fwl, ocont) swl, snorm = spt.sum_ec(fwl, onorm) - avgspecs.append(sp.replace('.spc.gz', '.rv.fits')) + avgspecs.append(sp.replace(".spc.gz", ".rv.fits")) nans, x = (np.isnan(scont), lambda z: z.nonzero()[0]) - sflx[nans] = 1. - scont[nans] = 1. + sflx[nans] = 1.0 + scont[nans] = 1.0 idx = np.where(scont <= 1e-3) - scont[idx] = 1. - sflx[idx] = 1. - final = sflx/scont + scont[idx] = 1.0 + sflx[idx] = 1.0 + final = sflx / scont idx = np.where((final > 0.986) & (final < 1.014)) pflx = final[idx] pwl = swl[idx] - pflx = np.interp(swl, pwl, savgol_filter(pflx, 3, 1, mode='nearest')) - final = final/pflx - spt.writeFits(final, swl*10, - savename=sp.replace('.spc.gz', '.rv.fits'), quiet=False, - externhd='../../{0}'.format(sp.replace('.spc.gz', '.fits')) ) + pflx = np.interp(swl, pwl, savgol_filter(pflx, 3, 1, mode="nearest")) + final = final / pflx + spt.writeFits( + final, + swl * 10, + savename=sp.replace(".spc.gz", ".rv.fits"), + quiet=False, + externhd="../../{0}".format(sp.replace(".spc.gz", ".fits")), + ) fig, ax = plt.subplots() ax.plot(swl, final) ax.set_ylim([0, 3]) - plt.savefig(figname=sp.replace('.spc.gz', '')) + plt.savefig(figname=sp.replace(".spc.gz", "")) plt.close(fig) -if os.path.exists('tmp.txt'): - os.system('rm tmp.txt') +if os.path.exists("tmp.txt"): + os.system("rm tmp.txt") spt.averagespecs(avgspecs) diff --git a/scripts/hdcompact.py b/scripts/hdcompact.py index 876343a..247ec33 100755 --- a/scripts/hdcompact.py +++ b/scripts/hdcompact.py @@ -18,25 +18,39 @@ __email__ = "dmfaes@gmail.com" -class MyParser(ArgumentParser): +class MyParser(ArgumentParser): def error(self, message): - sys.stderr.write('# ERROR! %s\n' % message) + sys.stderr.write("# ERROR! %s\n" % message) self.print_help() sys.exit(2) parser = MyParser(description=__doc__) -parser.add_argument('--version', action='version', - version='%(prog)s {0}'.format(__version__)) -parser.add_argument("-n", action="store", dest="n", - help=("Number of cores to be used [default: %(default)s]"), - type=int, default=2) - -group2 = parser.add_argument_group('other arguments') -group2.add_argument("-c", "--check", action="store_true", dest="chk", - help=("If this flag is enabled, it does not make new compressed files. It " +parser.add_argument( + "--version", action="version", version="%(prog)s {0}".format(__version__) +) +parser.add_argument( + "-n", + action="store", + dest="n", + help=("Number of cores to be used [default: %(default)s]"), + type=int, + default=2, +) + +group2 = parser.add_argument_group("other arguments") +group2.add_argument( + "-c", + "--check", + action="store_true", + dest="chk", + help=( + "If this flag is enabled, it does not make new compressed files. It " "will only removed the *.map* file if the " - "*.map.bz2 file exists"), default=False) + "*.map.bz2 file exists" + ), + default=False, +) args = parser.parse_args() @@ -44,23 +58,23 @@ def error(self, message): def compact(fp): if os.path.exists(fp): # print('# Creating '+fpnew) - fpnew = fp+'.bz2' + fpnew = fp + ".bz2" # fpnew = fp+'.gz' # with gzip.open(fpnew, 'wb') as fout: - with bz2.BZ2File(fpnew, 'wb', compresslevel=1) as fout: - fout.write(open(fp, 'rb').read()) + with bz2.BZ2File(fpnew, "wb", compresslevel=1) as fout: + fout.write(open(fp, "rb").read()) if os.path.exists(fpnew): if os.path.getsize(fpnew) > 40000: os.remove(fp) - print('# Processed '+fp) + print("# Processed " + fp) return 0 else: - print('# {0} was not found!'.format(fp)) + print("# {0} was not found!".format(fp)) return 1 -if __name__ == '__main__': - extensions = ['.tau', '.map'] +if __name__ == "__main__": + extensions = [".tau", ".map"] path = "." if len(sys.argv) == 2: @@ -75,18 +89,18 @@ def compact(fp): if not args.chk: clist.append(fp) else: - fpnew = fp+'.bz2' + fpnew = fp + ".bz2" if os.path.exists(fpnew): if os.path.getsize(fpnew) > 40000: os.remove(fp) - print('# Removing '+fp) + print("# Removing " + fp) if not args.chk and args.n > 1: pool = Pool(processes=args.n) # for result in pool.imap_unordered(compact, clist): result = pool.map(compact, clist) - print('# DONE!') + print("# DONE!") if not args.chk and args.n == 1: for c in clist: compact(c) - print('# DONE!') + print("# DONE!") diff --git a/scripts/norm_spec.py b/scripts/norm_spec.py index eef785c..37c1e7e 100755 --- a/scripts/norm_spec.py +++ b/scripts/norm_spec.py @@ -9,29 +9,33 @@ import numpy as np import os import sys -import pyfits as pf +import astropy.io.fits as pf import pyhdust.spectools as spt from pyhdust.tabulate import tabulate as tab from argparse import ArgumentParser -__version__ = "0.93" +__version__ = "0.94" __author__ = "Daniel Moser" __email__ = "dmfaes@gmail.com" -class MyParser(ArgumentParser): +class MyParser(ArgumentParser): def error(self, message): - sys.stderr.write('# ERROR! %s\n' % message) + sys.stderr.write("# ERROR! %s\n" % message) self.print_help() sys.exit(2) + parser = MyParser(description=__doc__) -parser.add_argument('--version', action='version', - version='%(prog)s {0}'.format(__version__)) -parser.add_argument("INPUT", help=("String to retrive the spectrum(a) " - "filename(s) (wildcards accepted)")) -# parser.add_argument("-l", "--line", action="store", dest="line", -# help=("Wavelength of the reference line [default: %(default)s]"), +parser.add_argument( + "--version", action="version", version="%(prog)s {0}".format(__version__) +) +parser.add_argument( + "INPUT", + help=("String to retrive the spectrum(a) " "filename(s) (wildcards accepted)"), +) +# parser.add_argument("-l", "--line", action="store", dest="line", +# help=("Wavelength of the reference line [default: %(default)s]"), # type=float, default=6562.79) args = parser.parse_args() @@ -40,11 +44,11 @@ def error(self, message): def norm_text(fname): f0 = np.loadtxt(fname) ny = spt.normalize_spec(f0[:, 0], f0[:, 1]) - outname = '{0}_norm{1}'.format(*os.path.splitext(fname)) + outname = "{0}_norm{1}".format(*os.path.splitext(fname)) # f0 = np.savetxt(outname, np.column_stack((f0[:, 0], ny))) - f1 = open(outname, 'w') - f1.writelines( tab(np.column_stack((f0[:, 0], ny)), tablefmt='plain') ) - print('# {0} saved!'.format(outname)) + f1 = open(outname, "w") + f1.writelines(tab(np.column_stack((f0[:, 0], ny)), tablefmt="plain")) + print("# {0} saved!".format(outname)) f1.close() return @@ -53,31 +57,33 @@ def norm_spec(fname): imfits = pf.open(fname) hdr = imfits[0].header flux = imfits[0].data - wl = np.arange(len(flux)) * hdr['CDELT1'] + hdr['CRVAL1'] + wl = np.arange(len(flux)) * hdr["CDELT1"] + hdr["CRVAL1"] nflux = spt.normalize_spec(wl, flux.astype(float)) - outname = '{0}_norm{1}'.format(*os.path.splitext(fname)) + outname = "{0}_norm{1}".format(*os.path.splitext(fname)) hdu = pf.PrimaryHDU(nflux) hdulist = pf.HDUList([hdu]) hdulist[0].header = hdr hdu.writeto(outname, overwrite=True) - print('# {0} saved!'.format(outname)) + print("# {0} saved!".format(outname)) imfits.close() return -if __name__ == '__main__': +if __name__ == "__main__": linput = args.INPUT for specfile in glob(linput): - if specfile.find('_norm.') > 0: - print('# Warning! Remove "_norm" from {0} if you want to normalize' - ' it!'.format(specfile)) + if specfile.find("_norm.") > 0: + print( + '# Warning! Remove "_norm" from {0} if you want to normalize' + " it!".format(specfile) + ) continue - if specfile.find('.fit') > 0: + if specfile.find(".fit") > 0: norm_spec(specfile) else: norm_text(specfile) if len(glob(linput)) == 0: - print('# ERROR! {0} files were not found!'.format(linput)) + print("# ERROR! {0} files were not found!".format(linput)) diff --git a/scripts/plot_spec.py b/scripts/plot_spec.py index a961fb2..4647fb9 100755 --- a/scripts/plot_spec.py +++ b/scripts/plot_spec.py @@ -12,68 +12,93 @@ import os import matplotlib.pyplot as plt -__version__ = "0.91" +__version__ = "0.92" __author__ = "Daniel Moser" __email__ = "dmfaes@gmail.com" -class MyParser(ArgumentParser): +class MyParser(ArgumentParser): def error(self, message): - sys.stderr.write('# ERROR! %s\n' % message) + sys.stderr.write("# ERROR! %s\n" % message) self.print_help() sys.exit(2) + parser = MyParser(description=__doc__) -parser.add_argument('--version', action='version', - version='%(prog)s {0}'.format(__version__)) -parser.add_argument("INPUT", help=("String to retrive the spectrum(a) " - "filename(s) (wildcards accepted wrapped with \"\")"), type=str) -parser.add_argument("-l", "--line", action="store", dest="line", - help=("Wavelength of the reference line [default: %(default)s]"), - type=float, default=6562.79) -parser.add_argument("-w", "--hwidth", action="store", dest="hwidth", - help=("Half-width of the line profile (in km/s) [default: %(default)s]"), - type=float, default=1000.) +parser.add_argument( + "--version", action="version", version="%(prog)s {0}".format(__version__) +) +parser.add_argument( + "INPUT", + help=( + "String to retrive the spectrum(a) " + 'filename(s) (wildcards accepted wrapped with "")' + ), + type=str, +) +parser.add_argument( + "-l", + "--line", + action="store", + dest="line", + help=("Wavelength of the reference line [default: %(default)s]"), + type=float, + default=6562.79, +) +parser.add_argument( + "-w", + "--hwidth", + action="store", + dest="hwidth", + help=("Half-width of the line profile (in km/s) [default: %(default)s]"), + type=float, + default=1000.0, +) # group2 = parser.add_argument_group('other arguments (overwrite power)') -# group2.add_argument("-r", "--remove", action="store_true", dest="rm", +# group2.add_argument("-r", "--remove", action="store_true", dest="rm", # help=("If this flag is enabled, it restores the spectrum(a) to the " # "original value (i.e., VELSHIFT = 0)"), default=False) -# group2.add_argument("-f", "--force-dv", action="store", dest="vs", +# group2.add_argument("-f", "--force-dv", action="store", dest="vs", # help=("Force VELSHIFT to the non-zero VS value (in km/s) " # "[default: None]"), type=float, default=0.) args = parser.parse_args() -if __name__ == '__main__': +if __name__ == "__main__": lbc = args.line linput = args.INPUT hwidth = args.hwidth - plev = 0. + plev = 0.0 lfiles = glob(linput) lfiles.sort() if len(lfiles): - fig, ax = plt.subplots(figsize=(8, 6*(1+len(lfiles)/25.))) + fig, ax = plt.subplots(figsize=(8, 6 * (1 + len(lfiles) / 25.0))) for fitsfile in lfiles: sdata = spt.loadfits(fitsfile) vl, fx = spt.lineProf(sdata[0], sdata[1], lbc=lbc, hwidth=hwidth) - ax.plot(vl, fx+plev) - ax.text(0.95*hwidth, 1.+plev+0.05/2, os.path.splitext(fitsfile)[0], - horizontalalignment='right', # transform=ax.transAxes, - verticalalignment='center', fontsize=8) + ax.plot(vl, fx + plev) + ax.text( + 0.95 * hwidth, + 1.0 + plev + 0.05 / 2, + os.path.splitext(fitsfile)[0], + horizontalalignment="right", # transform=ax.transAxes, + verticalalignment="center", + fontsize=8, + ) if plev == 0: ymin = min(fx) plev += 0.04 if len(glob(linput)) == 0: - print('# ERROR! {0} files were not found!'.format(linput)) + print("# ERROR! {0} files were not found!".format(linput)) else: - # ax.legend(loc='best', fancybox=True, framealpha=0.5, fontsize=8, + # ax.legend(loc='best', fancybox=True, framealpha=0.5, fontsize=8, # labelspacing=0.05) - ax.set_ylim([ymin, plev+max(fx)]) + ax.set_ylim([ymin, plev + max(fx)]) ax.set_xlim([-hwidth, hwidth]) - phc.savefig(fig, figname='plot_spec_'+phc.dtflag()) + phc.savefig(fig, figname="plot_spec_" + phc.dtflag()) diff --git a/scripts/shift_spec.py b/scripts/shift_spec.py index 58f3bfc..cf4db12 100755 --- a/scripts/shift_spec.py +++ b/scripts/shift_spec.py @@ -1,92 +1,127 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- -""" Program that shifts a FITS spectrum based on an line-profile Gaussian fit -(optimized to work within a few tenths of km/s). +"""Program that shifts a FITS spectrum based on an line-profile Gaussian +fit (optimized to work within a few tenths of km/s). """ -from glob import glob -import numpy as np -import sys -import pyhdust.spectools as spt -import pyhdust.phc as phc from argparse import ArgumentParser -import pyfits as pf -# from scipy.optimize import curve_fit +from glob import glob from matplotlib import pyplot as plt from warnings import warn -from lmfit import Model +import astropy.io.fits as pf +import numpy as np +import pyhdust.phc as phc +import pyhdust.spectools as spt +import sys + +try: + from lmfit import Model +except ImportError: + warn("lmfit module not installed!!!") -__version__ = "1.0" +__version__ = "1.1" __author__ = "Daniel Moser" __email__ = "dmfaes@gmail.com" -class MyParser(ArgumentParser): +class MyParser(ArgumentParser): def error(self, message): - sys.stderr.write('# ERROR! %s\n' % message) + sys.stderr.write("# ERROR! %s\n" % message) self.print_help() sys.exit(2) parser = MyParser(description=__doc__) -parser.add_argument('--version', action='version', - version='%(prog)s {0}'.format(__version__)) -parser.add_argument("INPUT", help=("String to retrive the spectrum(a) " - "filename(s) (wildcards accepted wrapped with \"\")"), type=str) -parser.add_argument("-l", "--line", action="store", dest="line", - help=("Wavelength of the reference line [default: %(default)s]"), - type=float, default=6562.79) -parser.add_argument("-w", "--half-width", action="store", dest="hw", - help=("Half-width in km/s of the reference line [default: %(default)s]"), - type=float, default=6562.79) - -group2 = parser.add_argument_group('other arguments (overwrite power)') -group2.add_argument("-r", "--remove", action="store_true", dest="rm", - help=("If this flag is enabled, it restores the spectrum(a) to the " - "original value (i.e., VELSHIFT = 0)"), default=False) -group2.add_argument("-f", "--force-dv", action="store", dest="vs", - help=("Force VELSHIFT to the non-zero VS value (in km/s) " - "[default: None]"), type=float, default=0.) +parser.add_argument( + "--version", action="version", version="%(prog)s {0}".format(__version__) +) +parser.add_argument( + "INPUT", + help=( + "String to retrive the spectrum(a) " + 'filename(s) (wildcards accepted wrapped with "")' + ), + type=str, +) +parser.add_argument( + "-l", + "--line", + action="store", + dest="line", + help=("Wavelength of the reference line [default: %(default)s]"), + type=float, + default=6562.79, +) +parser.add_argument( + "-w", + "--half-width", + action="store", + dest="hw", + help=("Half-width in km/s of the reference line [default: %(default)s]"), + type=float, + default=6562.79, +) + +group2 = parser.add_argument_group("other arguments (overwrite power)") +group2.add_argument( + "-r", + "--remove", + action="store_true", + dest="rm", + help=( + "If this flag is enabled, it restores the spectrum(a) to the " + "original value (i.e., VELSHIFT = 0)" + ), + default=False, +) +group2.add_argument( + "-f", + "--force-dv", + action="store", + dest="vs", + help=("Force VELSHIFT to the non-zero VS value (in km/s) " "[default: None]"), + type=float, + default=0.0, +) args = parser.parse_args() def gauss(x, a, x0, sigma): - y0=1. - return a*np.exp(-(x-x0)**2/(2*sigma**2))+y0 + y0 = 1.0 + return a * np.exp(-((x - x0) ** 2) / (2 * sigma**2)) + y0 def find_max_peak(vl, nflx, idx=-1): - """ Based on max flux. - """ + """Based on max flux.""" if idx == -1: - idx = len(vl)/2 + idx = len(vl) / 2 idx_ = phc.find_nearest(nflx[:idx], np.max(nflx[:idx]), idx=True) nfV, vlV = nflx[idx_], vl[idx_] idx_ = phc.find_nearest(nflx[idx:], np.max(nflx[idx:]), idx=True) - nfR, vlR = nflx[idx+idx_], vl[idx+idx_] + nfR, vlR = nflx[idx + idx_], vl[idx + idx_] return nfV, vlV, nfR, vlR def has_central_absorp(nf0, vlV, vlR, hw): - """ Return ``True`` or ``False`` if the line has a central absorption + """Return ``True`` or ``False`` if the line has a central absorption (i.e., it is an absorption line or it is an double peak emission line). """ if nf0 < 1: # print('# F(lb=0) < 1: A double-peak emission line was found!') return True - if vlR-vlV > hw/10: + if vlR - vlV > hw / 10: return True else: return False def gauss_fit(x, y, a0=None, x0=None, sig0=None, emission=True): - """ Return ``curve_fit``, i.e., ``popt, pcov``. + """Return ``curve_fit``, i.e., ``popt, pcov``. def gauss_fit(x, y, a0=None, x0=None, sig0=None, emission=True, ssize=0.05): - + # def gauss(x, a, x0, sigma): # y0=1. # return a*_np.exp(-(x-x0)**2/(2*sigma**2))+y0 @@ -138,98 +173,104 @@ def gauss_fit(x, y, a0=None, x0=None, sig0=None, emission=True, ssize=0.05): if x0 is None: x0 = x[np.where(y == func(y))] if sig0 is None: - sig0 = (np.max(x)-np.min(x))/10. + sig0 = (np.max(x) - np.min(x)) / 10.0 # if y0 is None: # y0 = np.median(y) gmodel = Model(gauss) - gmodel.set_param_hint('a', min=0.2, max=4) + gmodel.set_param_hint("a", min=0.2, max=4) if not emission: - gmodel.set_param_hint('a', min=-0.2, max=-4) - gmodel.set_param_hint('sigma', min=50, max=1000) + gmodel.set_param_hint("a", min=-0.2, max=-4) + gmodel.set_param_hint("sigma", min=50, max=1000) result = gmodel.fit(y, x=x, a=a0, x0=x0, sigma=sig0) fig, (ax0, ax1) = plt.subplots(2, 1) - ax0.plot(x, y, 'bo') - ax0.plot(x, result.init_fit, 'k--') - ax0.plot(x, result.best_fit, 'r-') + ax0.plot(x, y, "bo") + ax0.plot(x, result.init_fit, "k--") + ax0.plot(x, result.best_fit, "r-") ax0.set_title(fitsfile) print(lbc, np.min(x)) - idx = np.where((wl > lbc*0.98) & (wl < lbc*1.03)) - ax1.plot(wl[idx], flux[idx], 'o') + idx = np.where((wl > lbc * 0.98) & (wl < lbc * 1.03)) + ax1.plot(wl[idx], flux[idx], "o") plt.show(block=False) - cmd = phc.user_input('# Problem (y/other)? ') - if cmd.lower().startswith('y'): + cmd = phc.user_input("# Problem (y/other)? ") + if cmd.lower().startswith("y"): raise ValueError # phc.savefig(fig, figname=fitsfile) - return result.params['x0'] + return result.params["x0"] def apply_shift(imfits, vshift, wl=None): if wl is None: # print('watch!', imfits[0].header['CDELT1'], len(imfits[0].data)) - wl = np.arange(len(imfits[0].data)) * imfits[0].header['CDELT1'] + \ - imfits[0].header['CRVAL1'] + wl = ( + np.arange(len(imfits[0].data)) * imfits[0].header["CDELT1"] + + imfits[0].header["CRVAL1"] + ) - new_wl = wl*(1 + vshift/phc.c.cgs*1e5) + new_wl = wl * (1 + vshift / phc.c.cgs * 1e5) - imfits[0].header['CDELT1'] = (new_wl[-1]-new_wl[0])/(len(new_wl)-1) - imfits[0].header['CRVAL1'] = new_wl[0] + imfits[0].header["CDELT1"] = (new_wl[-1] - new_wl[0]) / (len(new_wl) - 1) + imfits[0].header["CRVAL1"] = new_wl[0] # imfits[0].data = np.interp(new_wl, wl, flux) old_sh = 0 - if 'VELSHIFT' in imfits[0].header: - old_sh = imfits[0].header['VELSHIFT'] - imfits[0].header['VELSHIFT'] = (vshift + old_sh, - 'Lambda shift in km/s') - - print('# {0} updated with {1:.1f} km/s shift (accum. {2:.0f} km/s)'. - format(fitsfile, vshift, vshift+old_sh)) + if "VELSHIFT" in imfits[0].header: + old_sh = imfits[0].header["VELSHIFT"] + imfits[0].header["VELSHIFT"] = (vshift + old_sh, "Lambda shift in km/s") + + print( + "# {0} updated with {1:.1f} km/s shift (accum. {2:.0f} km/s)".format( + fitsfile, vshift, vshift + old_sh + ) + ) return -if __name__ == '__main__': +if __name__ == "__main__": lbc = args.line linput = args.INPUT for fitsfile in glob(linput): - imfits = pf.open(fitsfile, mode='update') + imfits = pf.open(fitsfile, mode="update") if args.rm is True: - if 'VELSHIFT' in imfits[0].header: - vshift = -imfits[0].header['VELSHIFT'] + if "VELSHIFT" in imfits[0].header: + vshift = -imfits[0].header["VELSHIFT"] apply_shift(imfits, vshift) - else: - warn('No VELSHIFT for {0}'.format(fitsfile)) + else: + warn("No VELSHIFT for {0}".format(fitsfile)) elif args.vs != 0: v0 = 0 - if 'VELSHIFT' in imfits[0].header: - v0 = imfits[0].header['VELSHIFT'] - vshift = args.vs-v0 + if "VELSHIFT" in imfits[0].header: + v0 = imfits[0].header["VELSHIFT"] + vshift = args.vs - v0 apply_shift(imfits, vshift) else: try: flux = imfits[0].data - wl = np.arange(len(flux)) * imfits[0].header['CDELT1'] + \ - imfits[0].header['CRVAL1'] + wl = ( + np.arange(len(flux)) * imfits[0].header["CDELT1"] + + imfits[0].header["CRVAL1"] + ) vl, nflx = spt.lineProf(wl, flux, lbc=lbc, hwidth=args.hw) nfxsig = np.std(nflx) emission = True if np.percentile(nflx, 5) + nfxsig < 1: emission = False - if np.percentile(nflx, 95) - 1.5*nfxsig > 1: + if np.percentile(nflx, 95) - 1.5 * nfxsig > 1: emission = True vshift = -gauss_fit(vl, nflx, emission=emission) if np.abs(vshift) > args.hw: - raise ValueError('VelShift out of bounds!') + raise ValueError("VelShift out of bounds!") apply_shift(imfits, vshift, wl) except: - warn('Gaussian fit did not work for {0}'.format(fitsfile)) - vshift = phc.user_input('Enter a vshift value (km/s): ') + warn("Gaussian fit did not work for {0}".format(fitsfile)) + vshift = phc.user_input("Enter a vshift value (km/s): ") try: vshift = float(vshift) apply_shift(imfits, vshift, wl) @@ -239,4 +280,4 @@ def apply_shift(imfits, vshift, wl=None): imfits.close() if len(glob(linput)) == 0: - warn('{0} files were not found!'.format(linput)) + warn("{0} files were not found!".format(linput)) diff --git a/setup.py b/setup.py index 157e7c7..46495be 100755 --- a/setup.py +++ b/setup.py @@ -60,7 +60,7 @@ def recfiles(chdir, path): setup( name="pyhdust", - version="1.5.12", + version="1.5.12-1", description=( "Analysis tools for multi-technique astronomical data and hdust models" ),