[問題] 如何寫出系統中週期比和頻率比的關係,外加表現出波峰值

看板Fortran作者boto18 (暴投)時間15年前 (2009/04/30 14:20)推噓0(0推 0噓 18→)

留言18則, 3人參與討論串1/1

由於已經可以把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

→

YCTzeng

04/30 15:03, , 1^F

04/30 15:03, 1^F

→

charlesdc

04/30 18:15, , 2^F

04/30 18:15, 2^F

→

boto18

05/01 13:15, , 3^F

05/01 13:15, 3^F