| 1 |
#!/usr/bin/env python
|
| 2 |
"""convolver
|
| 3 |
|
| 4 |
Convolutes a set of computed frequencies (in cm^-1) with Lorentzian
|
| 5 |
linewidth to yield a simulated spectrum
|
| 6 |
|
| 7 |
Usage: convolver
|
| 8 |
|
| 9 |
Options:
|
| 10 |
-h, --help show this help
|
| 11 |
-f, --freqs=... use the specified frequencies (read from 2nd column)
|
| 12 |
-o, --out=... put the convoluted spectrum into this file
|
| 13 |
-w, --width=... use the specified Lorentzian linewidth (in cm^-1)
|
| 14 |
-n, --npoints=... use the specified number of data points
|
| 15 |
|
| 16 |
Example:
|
| 17 |
convolver -f FieldOff_Freq.data -o FieldOff_spectrum.dat -w 10
|
| 18 |
|
| 19 |
"""
|
| 20 |
|
| 21 |
__author__ = "Dan Gezelter (gezelter@nd.edu)"
|
| 22 |
__copyright__ = "Copyright (c) 2013 by the University of Notre Dame"
|
| 23 |
__license__ = "OpenMD"
|
| 24 |
|
| 25 |
import sys
|
| 26 |
import getopt
|
| 27 |
import string
|
| 28 |
import math
|
| 29 |
import random
|
| 30 |
from sets import *
|
| 31 |
import numpy
|
| 32 |
|
| 33 |
|
| 34 |
_haveFreqFileName = 0
|
| 35 |
_haveOutFileName = 0
|
| 36 |
_haveWidth = 0
|
| 37 |
_haveNpoints = 0
|
| 38 |
|
| 39 |
freqs = []
|
| 40 |
spectrum = []
|
| 41 |
|
| 42 |
def usage():
|
| 43 |
print __doc__
|
| 44 |
|
| 45 |
def lorentzian(x, background, amplitude, center, fwhm) :
|
| 46 |
# A lorentzian peak with:
|
| 47 |
# Constant Background : background
|
| 48 |
# Peak height above background : amplitude
|
| 49 |
# Central value : center
|
| 50 |
# Full Width at Half Maximum : fwhm
|
| 51 |
return background+(amplitude/numpy.pi)/(1.0+((x-center)/fwhm)**2)
|
| 52 |
|
| 53 |
|
| 54 |
def readFreqFile(FreqFileName):
|
| 55 |
FreqFile = open(FreqFileName, 'r')
|
| 56 |
|
| 57 |
for line in FreqFile:
|
| 58 |
L = line.split()
|
| 59 |
freq = float(L[1])
|
| 60 |
freqs.append(freq)
|
| 61 |
FreqFile.close()
|
| 62 |
freqsort = sorted(freqs)
|
| 63 |
min = freqsort[0]
|
| 64 |
max = freqsort[-1]
|
| 65 |
return (min, max)
|
| 66 |
|
| 67 |
def doSpectrum(min, max, width, nPoints):
|
| 68 |
df = float(max-min)/float(nPoints)
|
| 69 |
|
| 70 |
integral = 0.0
|
| 71 |
for i in range(nPoints):
|
| 72 |
f = min + df*i
|
| 73 |
spect = 0.0
|
| 74 |
for j in range(len(freqs)):
|
| 75 |
spect = spect + lorentzian(f, 0.0, 1.0, freqs[j], width)
|
| 76 |
integral += spect * df
|
| 77 |
spectrum.append( [f, spect] )
|
| 78 |
|
| 79 |
for i in range(len(spectrum)):
|
| 80 |
spectrum[i][1] = spectrum[i][1] / integral
|
| 81 |
|
| 82 |
|
| 83 |
def writeFile(OutFileName):
|
| 84 |
|
| 85 |
OutFile = open(OutFileName, 'w')
|
| 86 |
|
| 87 |
for i in range(len(spectrum)):
|
| 88 |
|
| 89 |
OutFile.write('%f\t%f\n' % (spectrum[i][0], spectrum[i][1]))
|
| 90 |
OutFile.close()
|
| 91 |
|
| 92 |
def main(argv):
|
| 93 |
try:
|
| 94 |
opts, args = getopt.getopt(argv, "hf:o:w:", ["help", "xyz=", "out=","width="])
|
| 95 |
except getopt.GetoptError:
|
| 96 |
usage()
|
| 97 |
sys.exit(2)
|
| 98 |
for opt, arg in opts:
|
| 99 |
if opt in ("-h", "--help"):
|
| 100 |
usage()
|
| 101 |
sys.exit()
|
| 102 |
elif opt in ("-f", "--freq"):
|
| 103 |
FreqFileName = arg
|
| 104 |
global _haveFreqFileName
|
| 105 |
_haveFreqFileName = 1
|
| 106 |
elif opt in ("-o", "--out"):
|
| 107 |
OutFileName = arg
|
| 108 |
global _haveOutFileName
|
| 109 |
_haveOutFileName = 1
|
| 110 |
elif opt in ("-w", "--width"):
|
| 111 |
width = float(arg)
|
| 112 |
global _haveWidth
|
| 113 |
_haveWidth = 1
|
| 114 |
elif opt in ("-n", "--npoints"):
|
| 115 |
nPoints = int(arg)
|
| 116 |
global _haveNpoints
|
| 117 |
_haveNpoints = 1
|
| 118 |
|
| 119 |
if (_haveFreqFileName != 1):
|
| 120 |
usage()
|
| 121 |
print "No frequency file was specified"
|
| 122 |
sys.exit()
|
| 123 |
|
| 124 |
if (_haveWidth != 1):
|
| 125 |
width = 1.0
|
| 126 |
if (_haveNpoints != 1):
|
| 127 |
nPoints = 1000
|
| 128 |
|
| 129 |
(min, max) = readFreqFile(FreqFileName)
|
| 130 |
print min, max
|
| 131 |
|
| 132 |
min = min - 5.0 * width
|
| 133 |
max = max + 5.0 * width
|
| 134 |
doSpectrum(min, max, width, nPoints)
|
| 135 |
|
| 136 |
if (_haveOutFileName == 1):
|
| 137 |
writeFile(OutFileName)
|
| 138 |
else:
|
| 139 |
usage()
|
| 140 |
print "No output file was specified"
|
| 141 |
sys.exit()
|
| 142 |
|
| 143 |
if __name__ == "__main__":
|
| 144 |
if len(sys.argv) == 1:
|
| 145 |
usage()
|
| 146 |
sys.exit()
|
| 147 |
main(sys.argv[1:])
|
