[Fortran] 纯文本查看 复制代码
!! Program to study the signal propagation in Hopf Media ( CSTR This File )
Program HBP_ICC_Main
Use MSIMSLMS
Use MSIMSLSS
Use MSFLIB
Use PORTLIB
Use Mymath_Module
Use Graph_Module
Use Model_Module
!! ------------------------------------------------------------------------------- !!
*4
Integer ji_Show( NS ) Integer, Parameter :: NS = 200, NV=3, M_Data = 50000---怎么构建循环，求300个，400个、、、1000个的？不需要多次重复输入。
Real x_Data( NS, M_Data ), x_Data1( M_Data ), x_PSD( 8193 )
Real x_Ary( NS, NV ), x_Ary0( NS, NV ), x_Gauss( NS )
Real xt(NV), dxdt(NV)
Integer ix, ix0, iy(NS), iy0(NS)
Real xx(NS)
Type (xyCOORD) xyPos
Character a1*3,a2*4 ,a3
Integer, Parameter :: N_C = 10, N_N = 20 !the number of coupling constants and noise
Real x_Matrix( N_C, N_N, NS, 4 ), xa_tmp( NS )
Real x_Result(7), xa_R( NS ), xa_Noise( N_N ), xa_CC( N_C )
!! Initialize the Screen for Graphic Drawing ------------------
Call InitGraphWindow( #FF0000 )
Call SetTextWindow( 2, 50, 80, 120 )
jfbk = SetTextColorRGB( 0 )
jfbk = SetColorRGB( #00FFFF )
!Call SetGraphWindow( 10, 10, 650, 510)
! For CSTR, T~2.0; So for DT=1E-3, T~2000*DT ; Set n_Tao=40, we have 50 points
! each period. For NT=2E6, we have 1000 periods. For 40000 data points, the
! total steps is 40000*40 = 1.6E6. So set n_Try=2E5, NT=2E6.
! For ICC, T~2.0; for DT=2E-3, n_Tao = 10, we have 100 points per period.
! Set NT=1E6, n_Try = 2E5, n_Data=40000, totally 400 periods ------
DT = 2E-3 ; NT = 6E5 ; n_tao = 10 ; n_try = 1E5
jf_ShowTime = 0 ; jf_Noise = 1 ; n_Delay = 1;jf_SaveTime = 0 !we can have a result via varing n_Delay
x_CP = .60 ; x_Max = 1.0 ; x_Min = 0.0
jf_BiDirection = 0 ! Birectional Flag !0(one_way),1(two_way)
n_Show = MIN( NS, 4 ) ; ji_Show(1:n_Show) = (/ 1, 2, 5, 10 /)!, 20, 30 /)
!! ------------------------------------------------------------------------------- !!
n_CC = 1 ; n_Noise = 11
! x_CC = 0.25 ; x_Noise = 0.04 ;
Do j_C = 1, n_CC ! For Coupling Constant
Do j_N = 1, n_Noise ! For Noise Intensity
!2.4
xa_CC( j_C ) = 3.0 + ( j_C-1 ) *2.0; xa_Noise( j_N ) = -2.0 + (j_N-1) * 0.2
x_CC = xa_CC( j_C ) ; x_Noise = 10.**( xa_Noise( j_N ) ) ! ~0.004 - 0.63
Call SetGraphWindow( 10, 10, 410, 810 )
j_t1 = Time() ! Start Time For Each Run -----
n_Step = 0 ; x_Gauss( 1:NS ) = 0.0 ; x_Ary = 0.1 ; n_Data = 0
!!----------------open a file to save the time series--------------------
if( jf_SaveTime == 1 ) then
write( a1, '(i3.3)' ) NS
write( a2, '(f4.1)' ) xa_Noise(j_N)
Write( a3, '(f4.1)' ) x_CC
Open( 5, file='Data\N_'//a1//'_noise_'//a2//'_CC_'//a3//'.dat', status='unknown' )
endif
!!---------------------------------------------------------
Do it = 1, NT
if( jf_Noise > 0 .AND. MOD( it-1, n_Delay ) == 0 ) then
Call RNNOA( jf_Noise, x_Gauss( 1:jf_Noise ) )
Do jj = 1, jf_Noise
if( x_Gauss(jj) < -5.0 ) x_Gauss(jj) = -5.0 ! To Avoid Overflow
if( x_Gauss(jj) > 5.0 ) x_Gauss(jj) = 5.0
Enddo
endif
x_Ary0 = x_Ary
Do js = 1, NS
! Using the ICC Model -------------------------------
xt(1:NV) = x_Ary( js, 1:NV )
xt_cp = x_CP * ( 1.+ x_Noise * x_Gauss(js) ) ! Mulplicative Noise to CP
Call Derivative( ICC_MODEL, it*DT, xt, dxdt, xt_CP )
if( js > 1 ) then
x_Ary(js,1) = x_Ary(js,1) + ( x_Ary0(js-1,1)-x_Ary0(js,1) ) * x_CC * DT ! Coupling
!x_Ary(js,1) = x_Ary(js,1) + x_Ary0(js-1,1 ) * x_CC * DT ! Coupling
endif
if( jf_BiDirection == 1 .AND. js < NS ) then
x_Ary(js,1) = x_Ary(js,1) + ( x_Ary0(js+1,1)-x_Ary0(js,1) ) * x_CC * DT ! Coupling
endif
!! Now Add the Effect of Noise ----------
!x_Ary( js,1 ) = x_Ary( js,1 ) + SQRT( 2.*x_Noise*DT ) * x_Gauss(js)
!x_Ary( js,1 ) = x_Ary( js,1 ) + x_Noise * x_Gauss(js) * DT ! Additive Noise
x_Ary( js, : ) = x_Ary( js, : ) + dxdt(:) * DT
Do j_V = 1, NV
if( x_Ary( js, j_V ) < 0 ) x_Ary( js, j_V ) = 0.0
if( x_Ary( js, j_V ) > 1 ) x_Ary( js, j_V ) = 1.0
Enddo
Enddo
if( it <= n_try ) cycle ! bypass the transient time
n_Step = n_Step + 1
if( n_Step == n_tao ) then
n_data = n_data + 1
x_Data( 1:NS, n_data ) = x_Ary(1:NS, 1)
!! -- Draw the graph on screen ---------------
if( jf_ShowTime == 1 ) then
ix = INT2(REAL(it-n_Try) / (NT-n_Try) * 410.)
Do ii = 1, n_Show
xx( ii ) = x_Ary( ji_Show(ii), 1 )
iy(ii) = 100*ii - INT2( REAL(xx(ii)-x_Min ) / (x_Max-x_Min) *100. )
Enddo
if( n_Data > 1 ) then
Do ii = 1, n_Show
Call Moveto( ix0, iy0(ii), xyPos )
jfbk = LineTo( ix, iy(ii) )
Enddo
endif
ix0 = ix ; iy0 = iy
endif
endif
if( n_Data >= M_Data ) goto 200
Enddo
200 Continue
! Now calculate the SNR Information and Save the Summary ------------------
Do js = 1, NS
x_Data1( 1:n_Data ) = x_Data( js, 1:n_Data )-----------主要是这里怎么求200个耦合细胞输出信噪比的平均值！
! 20 piece
Call Get_SNR_Data_2( n_Data, x_Data1, DT*n_Tao, 20, x_Result(1:7) )
x_Matrix( j_C, j_N, js, 1 ) = x_Result(1) ! Freq
x_Matrix( j_C, j_N, js, 2 ) = x_Result(2) ! Signal
x_Matrix( j_C, j_N, js, 3 ) = x_Result(3) ! Noise
x_Matrix( j_C, j_N, js, 4 ) = x_Result(6) ! SNR
xa_R( js ) = x_Result(6)
Enddo