| 1 |
jmarr |
4031 |
#!/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:])
|
