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:])
|