作者boto18 (暴投)
看板Fortran
標題[問題] 如何寫出系統中週期比和頻率比的關係,外加表現出波峰值
時間Thu Apr 30 14:20:01 2009
由於已經可以把data值跑出來但要如何取一段時間內
找出波峰的最大值, 外加求出週期比和頻率比的關係
程式如下fortran77語法有些真的很難懂
implicit real*8 (a-h,o-z)
parameter (maxnf=200000)
dimension f(maxnf),x(maxnf),v(maxnf),a(maxnf),res(maxnf)
open(file='2freq.out', unit=20)
PI=3.141592654
ncycle=100
F0=1.0
F1=0.5
phi1=0.0*PI
xi=0.0
ratio=2
w=1.0*2*PI
w2=w*ratio
period=(2.*PI)/w
dt=0.0025*period
nf=int(period*real(ncycle)/dt)
wn=2*PI
alpha=0.05
fy=1.0
Fs=F1
do 234 jj=1,nf
t=dt*real(jj-1)
f(jj)=F0*cos(w*t)+F1*cos(w2*t+phi1)
234 continue
call elastop(xi,wn,fy,alpha,nf,dt,f,x,v,res,xmax,
& resmax,xmu)
do 331 j=1,nf,1
t=dt*real(j-1)
tscale=t/period
acc=f(j)-2.*xi*wn*v(j)-res(j)
write(20,300) tscale,x(j),v(j),acc,res(j),f(j)
331 continue
300 format (8f9.4)
stop
end
subroutine elastop(xi,wn,fy,alpha,nf,dtorig,f,x,xd,res,xmax,
& resmax,xmu)
implicit real*8 (a-h,o-z)
dimension f(*),x(*),res(*),xd(*)
fyratio=0.005 ! 1% is allowed
toler=0.001 ! 0.001 in
vtol=wn*toler
maxiter=5
cp=1.0
xxmax=0.0
xresmax=0.0
xdd0=cp*f(1)
xkt1=(1.0-alpha)*wn*wn
xkt2=alpha*wn*wn
fy1=(1.0-alpha)*fy
xy=fy1/xkt1
x(1)=0.0
xd(1)=0.0
pcenter=0.0
nstate=0
xold=0.0
xdold=0.0
xddold=xdd0
b=2.0*xi*wn
t=0.0
do 100 i=1,nf
dtmove=0.0
tprev=dtorig*(i-1)
dt=dtorig
t=tprev+dtorig
df=f(i+1)-f(i)
fprev=f(i)
force=f(i+1)
goto 467
466 icut=1 ! reduce time step
dt=dt*0.5
df=df*0.5
force=fprev+df
t=tprev+dt
468 if(icut.eq.0) then
dt=dtorig-dtmove
t=tprev+dtorig
fprev=force
df=f(i+1)-fprev
force=f(i+1)
endif
467 continue
temp1=4.0/dt/dt + 2.0*b/dt
xtemp=xnew
do 234 iter=1,maxiter
if (nstate.eq.0) then
temp2=xkt1+xkt2
else
temp2=xkt2
endif
temp3=4.0/dt + b
temp5=force+temp1*xold +temp3*xdold+ xddold
if (nstate.eq.0) then
gtemp=(xtemp-pcenter)*xkt1+xtemp*xkt2
elseif (nstate.eq.1) then
gtemp=fy1+xtemp*xkt2
elseif (nstate.eq.2) then
gtemp=-fy1+xtemp*xkt2
endif
temp6=temp1*xtemp+gtemp
dx= (temp5-temp6)/(temp1+temp2)
xtemp=xtemp+dx
if (abs(dx) .gt.toler) then
goto 234
else
goto 355
endif
234 continue
355 continue
xnewtemp=xtemp
xddnewtemp = 4.0/dt/dt*(xnewtemp-xold) -4.0/dt*xdold -xddold
xdnewtemp = xdold +(xddold+xddnewtemp)*0.5*dt
subroutine linear2(xi,wn,nf,dt,f,x,v,a,xmax)
implicit real*8 (a-h,o-z)
dimension array1(2,2),array2(2,2)
dimension f(*),x(nf),v(nf),a(nf)
wd=wn*sqrt(1.0-xi**2)
dino1=exp((-1.0)*xi*wn*dt)
dino2=xi/sqrt(1.0-xi**2)
dino3=cos(wd*dt)
dino4=sin(wd*dt)
dino5=(2.0*xi**2-1.0)/(wn**2*dt)
dino6=(2.0*xi)/(wn**3*dt)
array1(1,1)=dino1*(dino3+dino2*dino4)
array1(1,2)=(dino1/wd)*dino4
array1(2,1)=((-1.0)*wn/sqrt(1.0-xi**2))*dino1*dino4
array1(2,2)=dino1*(dino3-dino2*dino4)
array2(1,1)=dino1*((dino5+xi/wn)*dino4/wd
& +(dino6+1./wn**2)*dino3)-dino6
array2(1,2)=(-1.)*dino1*(dino5*dino4/wd+dino6*dino3)
& -1./wn**2+dino6
array2(2,1)=(-1.)*(array1(1,1)-1.)/(wn**2*dt)-array1(1,2)
array2(2,2)=(-1.)*array2(2,1)-array1(1,2)
x(1)=0.0
v(1)=0.0
a(1)=f(1)
xold=x(1)
vold=v(1)
fold=f(1)
xmax=x(1)
do 10 i=2,nf
xnew=array1(1,1)*xold+array1(1,2)*vold
& -array2(1,1)*fold-array2(1,2)*f(i)
vnew=array1(2,1)*xold+array1(2,2)*vold
& -array2(2,1)*fold-array2(2,2)*f(i)
anew=f(i)-2.0*xi*wn*vnew-wn**2*xnew
x(i)=xnew
v(i)=vnew
a(i)=anew
fold=f(i)
xold=xnew
vold=vnew
if(abs(x(i)).gt.xmax) xmax=abs(x(i))
10 continue
return
end
subroutine whitegen(w,nw,dt,S0,seed)
implicit real*8 (a-h,o-z)
dimension w(*)
PI=3.14159
cons= sqrt(2*PI*S0/dt)
do 11 i=1,nw,2
call myrand(seed)
temp1=seed
call myrand(seed)
temp2=seed
w(i) =cons*sqrt(-2.0*dlog(temp1))*cos(2.0*PI*temp2)
w(i+1)=cons*sqrt(-2.0*dlog(temp1))*sin(2.0*PI*temp2)
11 continue
return
end
subroutine myrand (x)
implicit real*8(a-h,o-z)
m=67108864
rm=67108864.0
a=3125.0
temp=x*a*rm
xxrand=mod(temp,rm)
x=xxrand/rm
return
end
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.118.205.61
1F：→ YCTzeng:抱歉，根本看不懂你在問什麼 04/30 15:03
2F：→ charlesdc:建議先自己先整個看懂吧 才知道要編輯哪裡 04/30 18:15
3F：→ boto18:我再說清楚一些 05/01 13:15
4F：→ boto18:我已經跑出所需要的data值但現在需要作量化的圖形 05/01 13:16
5F：→ boto18:想要利用do loop的形式 取在特定時間內圖形中波峰最大值 05/01 13:18
6F：→ boto18:外加求出此最大波峰值對應的週期和相鄰的兩波峰的數值 05/01 13:19
7F：→ YCTzeng:特定時間是指某一次迴圈當中的意思嗎?還是你的data有時間? 05/01 21:13
8F：→ boto18:我的程式跑出來第一行就是時間,我現在只想取其中一個區段 05/02 15:27
9F：→ boto18:例如30~100之間,ncycle=100就是我所設定的最大值 05/02 15:29
10F：→ boto18:程式跑來的第二行是位移 我已經運用繪圖軟體畫出所需之圖形 05/02 15:30
11F：→ boto18:但由於點數過多無法從肉眼分辨圖形中的波峰最大值 05/02 15:31
12F：→ charlesdc:設一個矩陣去吃再利用矩陣函式去找最大值?? 05/02 17:20
13F：→ boto18:我之前就是用矩陣的情況去寫,結果教授不接受 05/02 17:51
14F：→ YCTzeng:多宣告一個變數,初值=0.0,進入迴圈後判斷波峰是否大於變數 05/02 18:38
15F：→ YCTzeng:,.true.的話,變數=波峰,等迴圈全部作完,就是波峰最大值, 05/02 18:40
16F：→ charlesdc:你有問為什麼不接受嗎..... 05/02 21:06
17F：→ boto18:Y大不好意思 由於圖形的最大值可能不只一個 05/02 21:39
18F：→ boto18:我要取的兩個最大的波峰之間的週期,還有相鄰波峰的週期 05/02 21:40