FreeSRS

 The shock response spectrum (SRS) allows to estimate transient acceleration signals in terms of a maximum response of a dynamic vibration system. Its origin probably dates from civil engineering and the problem of seismic excitations, i.e. earthquake response of buildings, see Clough Penzien 1975. The SRS is a spectrum, i.e. a curve giving the maximum acceleration response depending on frequency. It is derived from a 1-mass-spring-damper system being excited by a base excitation, i.e. an acceleration. It makes use of the nice fact, that the acceleration response of a one-mass-spring-damper system with base excitation does not depend on the mass, but only on eigenfrequency and damping. freesrs calculates the shock response spectrum (SRS) and is a small library with Python and Fortran files. The Fortran files consist of SmallwoodFilter.f and GenShockTimeHistories.f They are linked to Python by using the comfortable f2py tool. The algorithm used is from David O. Smallwood, see An improved recursive formula for calculating shock response spectra. Shock and Vibration Bulletin, No. 51, May 1981. For more background information you may also have a look on http://en.wikipedia.org/wiki/Shock_response_spectrum. The code was developed on a windows computer and compiled using Python(x,y) and f2py. There is also a pure C++ version which compiles with gcc. I intended to use the Fortran version together with Python, i.e. as the "fast" library. The program is very fast, all CPU-time-critical parts are written in Fortran 77. There are 2 modules (or libraries): 'SmallwoodFilter' is auto-generated with f2py (version:2). Functions: ypp = smallwoodfilter(model,fraction,xi,omega_n,ts,t,xpp,n) ypp = smallwoodfilterabsaccelmodelexact(xi,omega_n,ts,t,xpp,n) ypp = smallwoodfilterreldisplmodelexact(xi,omega_n,ts,t,xpp,n) yppm = recursiveequ(b0,b1,b2,a1s,a2s,xppm,xppm1,xppm2,yppm1,yppm2) b0,b1,b2,a1s,a2s = absaccelmodel(xi,omega_n,ts,t) b0,b1,b2,a1s,a2s = reldisplmodel(xi,omega_n,ts,t) ypp = smallwoodfilterapp(b0,b1,b2,a1s,a2s,xpp,n) b0,b1,b2,a1s,a2s = absaccelmodelapprox(xi,omega_n,ts) b0,b1,b2,a1s,a2s = reldisplmodelapprox(xi,omega_n,ts) 'GenShockTimeHistories' is auto-generated with f2py (version:2). Functions: t,a = gensquaredsinusshock(f_shock,fs,ampli,n) t,a = genrectangularshock(t_shock,fs,ampli,n) t = gentimevector(ts,n) noofconstlogfreq1 = noofconstlogfreq1(nfo,f1,f2) noofconstlogfreq2 = noofconstlogfreq2(nfo,f1,f2) f = compfreqconstlog(nfo,nf,f1,f2) Folder cpp It contains the cpp and exe files (debug and release version) compiled for windows. You may use the acc.inp file (half-sinus transient) for testing. I recommend using the f77 / Python version. There, all spectra in the time and the frequency domain are accessible to the user via the Python environment, e.g. when running ipython. This allows easy analysis of all kinds of SRS, i.e. positive, negative or primary and residual SRS. For generating the plots I used matplotlib (which is included in Python(x,y)). The code should also run under Linux or any other OS! Please, let me know when you compile the code for another OS that I can provide it to the community. Dependencies The following python packages are needed numpy PyLab Fortran compiler, e.g. gfortran, mingw32, ... There might be bugs! Please report bugs to dynamicsmec. Examples SRS of sinus shock:SRS of rectangular shock: Python code: #     Copyright (C) 2004-2013 by FreeDynamics #     http://sourceforge.net/projects/FreeSRS/ #     Licensed under the GNU GPL   from numpy import * from pylab import * import freesrs     print freesrs.SmallwoodFilter.__doc__ print freesrs.GenShockTimeHistories.__doc__   # half sine shock f_shock = 400. # Hz shock frequency fs = 50000. # Hz sampling frequency amp = 1.   t_shock = 1./f_shock/2.; ts = 1./fs; nl = int(t_shock/ts)+1 t,a = freesrs.GenShockTimeHistories.gensquaredsinusshock(f_shock,fs,amp,nl) figure(1) plot(t,a) xlabel('sec') ylabel('m/s**2') f1 = 10.; f2 = 4000.; nfo = 30 nf = freesrs.GenShockTimeHistories.noofconstlogfreq1(nfo,f1,f2) #print nf frev = freesrs.GenShockTimeHistories.compfreqconstlog(nfo,nf,f1,f2) #print frev   xi = 0.01 srs = arange(float(len(frev))) for i in range(len(frev)):     omega_n = 2.*pi*frev[i]     ypp = freesrs.SmallwoodFilter.smallwoodfilter('a',0.1,xi,omega_n,ts,t,a,nl)     max_1 = max(ypp)     # determine ypp(tend) and yp(tend) if x0pp = 0 by num integration     app = ypp[nl-1]; ap = 0.5*(ypp[nl-1]+ypp[nl-2])*ts     max_2 = freesrs.AnalysePostVib(ap, app, xi, omega_n)     #print frev[i], max_1, max_2     if max_1 > max_2:         srs[i] = max_1         #print frev[i]     else:         srs[i] = max_2         #print frev[i], max_1, max_2   figure(2) plot (frev, srs) xlabel('Hz') ylabel('m/s**2') show()