1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85

	    import numpy as np
	    import math
	    import os,glob
	    import string
	    import matplotlib.pyplot as plt
	    from pylab import *
	    import datetime
	    starttime = datetime.datetime.now()
	    def computepga(Accelerate,LN,damp,Freq,nFreq):
	                  n =LN
	                  dt = 0.02
	                  zeta = damp
	                  nfreq = nFreq
	                  freq = Freq
	                  MAcc = [0]*nfreq 
	                  absacc = [0]*n
	                  for s in range(n):
	                                absacc[s] = abs( string.atof(Accelerate[s]) ) 
	                  for k in range(nfreq):
	                                Damfreq = freq[k] * sqrt( 1-zeta**2 )
	                                e_t  = math.exp( - zeta * dt * freq[k] )
	                                ss = math.sin( Damfreq * dt )
	                                cc = math.cos( Damfreq * dt )
	                                d_f = ( 2 * zeta**2 - 1 ) / ( (freq[k])**2 * dt )
	                                d_3t = zeta / ( freq[k]**3 * dt )
	                                A00= e_t*( ss * zeta/sqrt(1-zeta*zeta)+cc )
	                                A01= e_t*( ss/Damfreq )
	                                A10 =  -freq[k] * e_t * ss/sqrt(1-zeta*zeta)
	                                A11 =   e_t * ( -ss * zeta / sqrt(1-zeta*zeta) + cc )
	                                B00 =   e_t * (( d_f+zeta / freq[k] ) * ss/Damfreq+( 2*d_3t+1 / freq[k]**2) * cc )-2*d_3t
	                                B01 =  -e_t * ( d_f * ss/Damfreq + 2*d_3t*cc ) - 1/freq[k]**2 + 2*d_3t
	                                B10  =  e_t * ( (d_f +zeta/freq[k] ) * (cc - zeta/sqrt(1-zeta**2) *ss ) - ( 2*d_3t+1/freq[k]**2 ) * ( Damfreq*ss+zeta*freq[k]*cc) )+1/( freq[k]**2*dt )
	                                B11 =   e_t * ( 1/(freq[k]**2*dt )*cc+ss*zeta/( freq[k]*Damfreq*dt ))-1/( freq[k]**2*dt )           
	                                Displace = [0]*(n+1)
	                                Velocity = [0]*(n+1)
	                                AbsAcce = [0]*(n+1)
	                                absabsAcc = [0]*(n+1)
	                                for j in range(n-1):
	                                              Displace[j+1] = A00*Displace[j]+A01*Velocity[j]+B11*float(Accelerate[j])+B01*float(Accelerate[j+1]) 
	                                              Velocity[j+1] = A10*Displace[j]+A11*Velocity[j]+B10*float(Accelerate[j])+B11*float(Accelerate[j+1]) 
	                                              AbsAcce[j+1]  = -2*zeta*freq[k]*Velocity[j+1] -freq[k]**2*Displace[j+1]
	                                              absabsAcc[j+1] = abs(AbsAcce[j+1])
	                                if freq[k]>=10**10:
	                                              MAcc[k] = max(absacc)
	                                else:
	                                              MAcc[k] = max(absabsAcc)
	                  return MAcc
	    accfile = glob.glob('xj30bar/ACC/*S001*.acc')
	    numsite = len(accfile)
	    print "the number of site is: ",numsite
	    damp = 0.05
	    Freq1 = 0.1
	    Freq2 = 50
	    nFreq = 100
	    Frinc = math.log(Freq2/Freq1)/(nFreq-1)
	    Freq = [0]*nFreq
	    T = [0]*nFreq
	    Freq[0] = Freq1
	    for k in range(nFreq):
	                 Freq[k] = Freq1*math.exp(k*Frinc)
	                 T[k] = 1/Freq[k]
	    PGA = [[0]*nFreq]*numsite
	    for idx in range(numsite):
	                filename = accfile[idx]
	                file =  open(filename,'r')
	                linesList = file.readlines()
	                linesList = [line.strip().split() for line in linesList]
	                file.close()
	                tm = [x[0] for x in linesList[17:65552]]
	                acc = [x[1] for x in linesList[17:65552]]
	                LN =len(acc)
	                print "computing", filename
	                PGA = computepga(acc,LN,damp,Freq,nFreq)
	                plt.figure(1)       
	                plt.loglog(T,PGA,label = '%d' %(idx+1))
	                plt.grid()
	                plt.title("PGA")
	                plt.xlabel("T")
	                plt.ylabel("PGA")
	                plt.legend()
	    print 'plot figure now'     
	    endtime = datetime.datetime.now()
	    print (endtime - starttime).seconds,'seconds'
	    plt.show()
	    print('happy end')