API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jun 17th, 2017
82
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 24.09 KB | None | 0 0
raw download clone embed print report
  1. ! Chuong trinh se tinh n(t), p(t) và u(t)
  2. ! tham so dau vao: e(k), f0   (cac tham so nay se dc nhap vao, hoac viet mot truong chinh tinh toan rieng.
  3. ! Cac gia tri ban dau: n(t=0)=0 ; p(t=0) = 0; u(t=0) =0; t0=0;
  4.  
  5. Module global      ! Khai bao cac hang so se su dung trong chuong trinh
  6.  implicit none
  7.  save
  8.  real, parameter :: pi= 3.141593
  9.  real, parameter :: hbar = 658.d-3
  10.  complex, parameter :: iu = (0.,1.) , ui = (1.,0.)
  11.  real(8) :: t0
  12.  integer, parameter, public:: estep = 10  ! Luoi' k (cac gia tri cua k)
  13.  complex, dimension(2,2) :: sigma_z, sigma_y  ! khai bao 2 ma tran sigma
  14.  
  15. ! Pho gia tri nang luong theo k.
  16.  real(8) ::   f0, h
  17. real(8),  dimension(estep) ::e
  18. End module global
  19.  
  20. !============================================================================================================
  21.  
  22.  
  23. Module Runge                                                ! Runge-kutta gan dung bac 4.
  24.  use global
  25.  implicit none
  26.  contains                                                   ! Lenh nay dung de bao lai cac subroutine con cua subroutine Runge
  27.    subroutine rk4(t, h, e,f0, p, n, u,A_in, B_in, C_in, D_in, E_in, F_in,G_in, H_in,I_in, J_in, K_in, L_in,&
  28.                               A_out, B_out, C_out, D_out, E_out, F_out, G_out, H_out, I_out, J_out, K_out, L_out, &
  29.                                t_out, p_out, n_out, u_out)
  30.                      
  31.  
  32.     implicit none
  33.     real(8), intent(in), dimension(estep) ::e
  34.     complex(8), intent(in) :: p                                ! Vi p là ham chi phu thuoc t, ko phu thuoc k.
  35.     complex(8), intent(out):: p_out
  36.  
  37.     complex(8), intent(in) , dimension(estep) :: n             ! mang n gom estep phan tu, moi phan tu ung voi mot gia tri cua k.
  38.     complex(8), intent(out), dimension(estep) :: n_out
  39.  
  40.     real(8), intent(in) :: t ,  h  , f0
  41.     real(8), intent(out) :: t_out                        ! t la thoi gian, h la buoc nhay cua t
  42.    
  43.     complex(8), intent (in), dimension(2,2,estep):: u
  44.     complex(8), intent (out), dimension(2,2,estep):: u_out  !Ung voi moi gia tri cua k, ham U la mot matran 2x2
  45.  
  46.     complex(8), intent(in),dimension(2,2,2,estep,estep,estep) :: A_in, B_in, C_in, D_in, E_in, F_in,&
  47.                                                                   G_in, H_in, I_in, J_in, K_in, L_in
  48.     complex(8), intent(out),dimension(2,2,2,estep, estep, estep) ::A_out, B_out, C_out, D_out, E_out, F_out, G_out, H_out, &
  49.                                                  I_out, J_out, K_out, L_out
  50.  
  51.  
  52.     ! Khai bao cac he so tam thoi trong subroutine:
  53.     real(8) :: t1, t3
  54.     complex(8) :: kp1, kp2, kp3, kp4, p1, p2, p3
  55.     complex(8), dimension(estep) :: kn1, kn2, kn3, kn4, n1, n2, n3
  56.     complex(8), dimension(2,2,estep):: ku1, ku2, ku3, ku4, u1, u2, u3
  57.     complex(8), dimension(2,2,2, estep, estep, estep) :: A_temp1, A_temp2, A_temp3, A_temp4, B_temp1, B_temp2, B_temp3, B_temp4,  &
  58.                                                          C_temp1, C_temp2, C_temp3, C_temp4, D_temp1, D_temp2, D_temp3, D_temp4  
  59.  
  60.  
  61.     complex(8), dimension(2,2,2, estep, estep, estep) ::  E_temp1, E_temp2, E_temp3, E_temp4, F_temp1, F_temp2, F_temp3, F_temp4, &
  62.                                                    G_temp1, G_temp2, G_temp3, G_temp4, H_temp1, H_temp2, H_temp3, H_temp4, &
  63.                                                    I_temp1, I_temp2, I_temp3, I_temp4, J_temp1, J_temp2, J_temp3, J_temp4, &
  64.                                                    K_temp1, K_temp2, K_temp3, K_temp4, L_temp1, L_temp2, L_temp3, L_temp4
  65.    
  66.     complex(8), dimension(2,2,2,estep,estep,estep):: A_1, A_2, A_3, B_1, B_2, B_3, C_1, C_2, C_3, D_1, D_2, D_3, &
  67.                                                      E_1, E_2, E_3, F_1, F_2, F_3, G_1, G_2, G_3, H_1, H_2, H_3, &
  68.                                                      I_1, I_2, I_3, J_1, J_2, J_3, K_1, K_2, K_3, L_1, L_2, L_3                                                      
  69.              
  70.  
  71.     call u_deriv(t,e,p,n,u,ku1)
  72.  
  73.     Call A_deriv(t,p,n,u,A_temp1)
  74.     call B_deriv(t,p,n,u,B_temp1)
  75.     call C_deriv(t,p,n,u,C_temp1)
  76.     call D_deriv(t,p,n,u,D_temp1)
  77.  
  78.     call p_deriv(t,p,n,u, A_in, B_in, C_in, D_in,kp1)
  79.  
  80.     call E_deriv(t,p,n,u,E_temp1)
  81.     call F_deriv(t,p,n,u,F_temp1)
  82.     call G_deriv(t,p,n,u,G_temp1)
  83.     call H_deriv(t,p,n,u,H_temp1)
  84.     call I_deriv(t,p,n,u,I_temp1)
  85.     call J_deriv(t,p,n,u,J_temp1)
  86.     call K_deriv(t,p,n,u,K_temp1)
  87.     call L_deriv(t,p,n,u,L_temp1)
  88.  
  89.     call n_deriv(t,p,n,u, E_1, F_1, G_1, H_1, I_1, J_1, K_1, L_1, kn1)
  90.  
  91.       t1  = t+ h/2.0
  92.       u1 = u + h*ku1/2.0
  93.       n1 = n + h*kn1/2.0
  94.       p1 = p + h*kp1/2.0
  95.  
  96.       A_1 = A_in + h*A_temp1/2.0
  97.       B_1 = B_in + h*B_temp1/2.0
  98.       C_1 = C_in + h*C_temp1/2.0
  99.       D_1 = D_in + h*D_temp1/2.0
  100.       E_1 = E_in + h*E_temp1/2.0
  101.       F_1 = F_in + h*F_temp1/2.0
  102.       G_1 = G_in + h*G_temp1/2.0
  103.       H_1 = H_in + h*H_temp1/2.0
  104.       I_1 = I_in + h*I_temp1/2.0
  105.       J_1 = J_in + h*J_temp1/2.0
  106.       K_1 = K_in + h*K_temp1/2.0
  107.       L_1 = L_in + h*L_temp1/2.0
  108.  
  109.  
  110.  
  111.     call u_deriv(t1, e, p1, n1, u1, ku2)
  112.  
  113.     Call A_deriv(t1,p1,n1,u1,A_temp2)
  114.     call B_deriv(t1,p1,n1,u1,B_temp2)
  115.     call C_deriv(t1,p1,n1,u1,C_temp2)
  116.     call D_deriv(t1,p1,n1,u1,D_temp2)
  117.  
  118.     call p_deriv(t1,p1,n1,u1, A_1, B_1, C_1, D_1,kp2)
  119.  
  120.     call E_deriv(t1,p1,n1,u1,E_temp2)
  121.     call F_deriv(t1,p1,n1,u1,F_temp2)
  122.     call G_deriv(t1,p1,n1,u1,G_temp2)
  123.     call H_deriv(t1,p1,n1,u1,H_temp2)
  124.     call I_deriv(t1,p1,n1,u1,I_temp2)
  125.     call J_deriv(t1,p1,n1,u1,J_temp2)
  126.     call K_deriv(t1,p1,n1,u1,K_temp2)
  127.     call L_deriv(t1,p1,n1,u1,L_temp2)
  128.  
  129.     call n_deriv(t1,p1,n1,u1, E_2, F_2, G_2, H_2, I_2, J_2, K_2, L_2, kn2)
  130.      
  131.       u2 = u + h*ku2/2.0
  132.       n2 = n + h*kn2/2.0
  133.       p2 = p + h*kp2/2.0
  134.  
  135.       A_2 = A_in + h*A_temp2/2.0
  136.       B_2 = B_in + h*B_temp2/2.0
  137.       C_2 = C_in + h*C_temp2/2.0
  138.       D_2 = D_in + h*D_temp2/2.0
  139.       E_2 = E_in + h*E_temp2/2.0
  140.       F_2 = F_in + h*F_temp2/2.0
  141.       G_2 = G_in + h*G_temp2/2.0
  142.       H_2 = H_in + h*H_temp2/2.0
  143.       I_2 = I_in + h*I_temp2/2.0
  144.       J_2 = J_in + h*J_temp2/2.0
  145.       K_2 = K_in + h*K_temp2/2.0
  146.       L_2 = L_in + h*L_temp2/2.0     
  147.      
  148.  
  149.  
  150.     call u_deriv(t1, e, p2, n2, u2, ku3)
  151.  
  152.     Call A_deriv(t1,p2,n2,u2,A_temp3)
  153.     call B_deriv(t1,p2,n2,u2,B_temp3)
  154.     call C_deriv(t1,p2,n2,u2,C_temp3)
  155.     call D_deriv(t1,p2,n2,u2,D_temp3)
  156.  
  157.     call p_deriv(t1, p2, n2, u2, A_2, B_2, C_2, D_2, kp3)
  158.  
  159.     call E_deriv(t1,p2,n2,u2,E_temp3)
  160.     call F_deriv(t1,p2,n2,u2,F_temp3)
  161.     call G_deriv(t1,p2,n2,u2,G_temp3)
  162.     call H_deriv(t1,p2,n2,u2,H_temp3)
  163.     call I_deriv(t1,p2,n2,u2,I_temp3)
  164.     call J_deriv(t1,p2,n2,u2,J_temp3)
  165.     call K_deriv(t1,p2,n2,u2,K_temp3)
  166.     call L_deriv(t1,p2,n2,u2,L_temp3)
  167.  
  168.     call n_deriv(t1, p2, n2, u2, E_2, F_2, G_2, H_2, I_2, J_2, K_2, L_2, kn3)
  169.  
  170.      t3 = t + h
  171.      u3 = u + h*ku3
  172.      n3 = n + h*kn3
  173.      p3 = p + h*kp3
  174.      A_3 = A_in + h*A_temp3
  175.      B_3 = B_in + h*B_temp3
  176.      C_3 = C_in + h*C_temp3
  177.      D_3 = D_in + h*D_temp3
  178.      E_3 = E_in + h*E_temp3
  179.      F_3 = F_in + h*F_temp3
  180.      G_3 = G_in + h*G_temp3
  181.      H_3 = H_in + h*H_temp3
  182.      I_3 = I_in + h*I_temp3
  183.      J_3 = J_in + h*J_temp3
  184.      K_3 = K_in + h*K_temp3
  185.      L_3 = L_in + h*L_temp3
  186.  
  187.      
  188.  
  189.     call u_deriv(t3, e, p3, n3, u3, ku4)
  190.  
  191.     Call A_deriv(t3,p3,n3,u3,A_temp4)
  192.     call B_deriv(t3,p3,n3,u3,B_temp4)
  193.     call C_deriv(t3,p3,n3,u3,C_temp4)
  194.     call D_deriv(t3,p3,n3,u3,D_temp4)
  195.  
  196.     call p_deriv(t3, p3, n3, u3, A_3, B_3, C_3, D_3, kp4)
  197.  
  198.     call E_deriv(t3,p3,n3,u3,E_temp4)
  199.     call F_deriv(t3,p3,n3,u3,F_temp4)
  200.     call G_deriv(t3,p3,n3,u3,G_temp4)
  201.     call H_deriv(t3,p3,n3,u3,H_temp4)
  202.     call I_deriv(t3,p3,n3,u3,I_temp4)
  203.     call J_deriv(t3,p3,n3,u3,J_temp4)
  204.     call K_deriv(t3,p3,n3,u3,K_temp4)
  205.     call L_deriv(t3,p3,n3,u3,L_temp4)
  206.  
  207.     call n_deriv(t3, p3, n3, u3, E_3, F_3, G_3, H_3, I_3, J_3, K_3, L_3, kn4)
  208.  
  209.      t_out = t3
  210.      u_out = u + h*(ku1 +ku4 + 2.0*(ku2 + ku3))/6.0
  211.      n_out = n + h*(kn1 +kn4 + 2.0*(kn2 + kn3))/6.0
  212.      p_out = p + h*(kp1 +kp4 + 2.0*(kp2 + kp3))/6.0
  213.  
  214.      A_out = A_in + h*(A_temp1 + A_temp4 +2.0*(A_temp2 + A_temp3))/6.0
  215.      B_out = B_in + h*(B_temp1 + B_temp4 +2.0*(B_temp2 + B_temp3))/6.0
  216.      C_out = C_in + h*(C_temp1 + C_temp4 +2.0*(C_temp2 + C_temp3))/6.0
  217.      D_out = D_in + h*(D_temp1 + D_temp4 +2.0*(D_temp2 + D_temp3))/6.0
  218.      E_out = E_in + h*(E_temp1 + E_temp4 +2.0*(E_temp2 + E_temp3))/6.0
  219.      F_out = F_in + h*(F_temp1 + F_temp4 +2.0*(F_temp2 + F_temp3))/6.0
  220.      G_out = G_in + h*(G_temp1 + G_temp4 +2.0*(G_temp2 + G_temp3))/6.0
  221.      H_out = H_in + h*(H_temp1 + H_temp4 +2.0*(H_temp2 + H_temp3))/6.0
  222.      I_out = I_in + h*(I_temp1 + I_temp4 +2.0*(I_temp2 + I_temp3))/6.0
  223.      J_out = J_in + h*(J_temp1 + J_temp4 +2.0*(J_temp2 + J_temp3))/6.0
  224.      K_out = K_in + h*(K_temp1 + K_temp4 +2.0*(K_temp2 + K_temp3))/6.0
  225.      L_out = L_in + h*(L_temp1 + L_temp4 +2.0*(L_temp2 + L_temp3))/6.0
  226.  
  227.     end subroutine rk4
  228.  
  229.  
  230.    Subroutine u_deriv(t,e,p,n,u,ku)
  231.  
  232.     implicit none
  233.  
  234.     real(8), intent(in) :: t
  235.     complex(8), intent(in) ::p
  236.     complex(8), intent(in), dimension(estep) :: n
  237.     complex(8), intent(in), dimension(2,2,estep) :: u
  238.     complex(8), intent(out), dimension(2,2,estep) ::ku
  239.     real(8), intent(in), dimension(estep) ::e
  240.  
  241.     integer :: i
  242.  
  243.     complex(8), dimension(2,2, estep) :: H                         !H la hamiltonian
  244.     complex(8), dimension(2,2) ::H_temp , u_temp , ku_temp    
  245.                              
  246.      Do i =1 , estep                                             ! Thiet lap H, moi phan tu cua H theo k la mot matran 2x2
  247.       H(:,:,i) = (e(i) + f0*n(i))* sigma_z + iu*f0*p*p*sigma_y
  248.      end do
  249.  
  250.      Do i = 1, estep
  251.  
  252.       ku(:,:,i) = -iu* matmul(H(:,:,i),u(:,:,i))                 !Day la pt A.11 , Lenh matmul() la nhan ma tran
  253.      end do
  254.  
  255.     end subroutine u_deriv
  256.  
  257.  
  258.     Subroutine A_deriv(t, p, n, u, A_temp)
  259.      implicit none
  260.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: A_temp
  261.       real(8), intent(in) :: t
  262.       complex(8), intent(in):: p
  263.       complex(8), dimension(estep), intent(in) :: n
  264.       complex(8), intent(in), dimension(2,2,estep) :: u
  265.       integer :: b,d,m,k,l,q
  266.  
  267.        Do b=1, 2
  268.         Do d=1, 2
  269.          Do m=1, 2
  270.           Do k=1, estep
  271.            Do q=1, estep
  272.             l = k - q
  273.             If (l >= 1) then
  274.             A_temp(b,d,m,k,l,q) = conjg(u(b,1,k))*u(d,1,l)*u(m,1,q)*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*conjg(p)
  275.             endif
  276.            end do
  277.           end do
  278.          end do
  279.         end do
  280.        end do
  281.      End subroutine A_deriv
  282.  
  283.  
  284.  
  285.      Subroutine B_deriv(t, p, n, u, B_temp)
  286.       implicit none
  287.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: B_temp
  288.       real(8), intent(in) :: t
  289.       complex(8), intent(in):: p
  290.       complex(8), dimension(estep), intent(in) :: n
  291.       complex(8), intent(in), dimension(2,2,estep) :: u
  292.       integer :: b,d,m,k,l,q
  293.        
  294.        Do b=1, 2
  295.         Do d=1, 2
  296.          Do m=1, 2
  297.           Do k=1, estep
  298.            Do q=1, estep
  299.             l = k + q
  300.             If (l >= 1 .and. l <= estep) then
  301.             B_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l) + n(l)*n(q)+n(l)*n(k) - n(k)*n(q))*conjg(p)
  302.             endif
  303.            end do
  304.           end do
  305.          end do
  306.         end do
  307.        end do
  308.      End subroutine B_deriv
  309.  
  310.  
  311.     Subroutine C_deriv(t, p, n, u, C_temp)
  312.      implicit none
  313.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: C_temp
  314.       real(8), intent(in) :: t
  315.       complex(8), intent(in):: p
  316.       complex(8), dimension(estep), intent(in) :: n
  317.       complex(8), intent(in), dimension(2,2,estep) :: u
  318.       integer :: b,d,m,k,l,q
  319.  
  320.        Do b=1, 2
  321.         Do d=1, 2
  322.          Do m=1, 2
  323.           Do k=1, estep
  324.            Do q=1, estep
  325.             l = k - q
  326.             If (l >= 1) then
  327.             C_temp(b,d,m,k,l,q) = conjg(u(b,1,k))*u(d,1,l)*u(m,1,q)*(n(k)+n(k)*n(l)+n(k)*n(q)-n(l)*n(q))*p
  328.             endif
  329.            end do
  330.           end do
  331.          end do
  332.         end do
  333.        end do
  334.      End subroutine C_deriv
  335.  
  336.  
  337.  
  338.     Subroutine D_deriv(t, p, n, u, D_temp)
  339.      implicit none
  340.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: D_temp
  341.       real(8), intent(in) :: t
  342.       complex(8), intent(in):: p
  343.       complex(8), dimension(estep), intent(in) :: n
  344.       complex(8), intent(in), dimension(2,2,estep) :: u
  345.       integer :: b,d,m,k,l,q
  346.  
  347.        Do b=1, 2
  348.         Do d=1, 2
  349.          Do m=1, 2
  350.           Do k=1, estep
  351.            Do q=1, estep
  352.             l = k + q
  353.             If (l >= 1 .and. l <= estep) then
  354.             D_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*p
  355.             endif
  356.            end do
  357.           end do
  358.          end do
  359.         end do
  360.        end do
  361.      End subroutine D_deriv
  362.  
  363.  
  364.  
  365.     Subroutine E_deriv(t, p, n, u, E_temp)
  366.      implicit none
  367.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: E_temp
  368.       real(8), intent(in) :: t
  369.       complex(8), intent(in):: p
  370.       complex(8), dimension(estep), intent(in) :: n
  371.       complex(8), intent(in), dimension(2,2,estep) :: u
  372.       integer :: b,d,m,k,l,q
  373.  
  374.        Do b=1, 2
  375.         Do d=1, 2
  376.          Do m=1, 2
  377.           Do k=1, estep
  378.            Do q=1, estep
  379.             l = k - q
  380.             If (l >= 1) then
  381.             E_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*conjg(u(m,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*p
  382.             endif
  383.            end do
  384.           end do
  385.          end do
  386.         end do
  387.        end do
  388.      End subroutine E_deriv
  389.  
  390.  
  391. !thuc ra F_deriv chi la conjg cua E_deriv, co the luoc bo cac giai doan tinh cac F,H,J,L nay de lam chuong trinh nhe hon.
  392.     Subroutine F_deriv(t, p, n, u, F_temp)
  393.      implicit none
  394.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: F_temp
  395.       real(8), intent(in) :: t
  396.       complex(8), intent(in):: p
  397.       complex(8), dimension(estep), intent(in) :: n
  398.       complex(8), intent(in), dimension(2,2,estep) :: u
  399.       integer :: b,d,m,k,l,q
  400.  
  401.        Do b=1, 2
  402.         Do d=1, 2
  403.          Do m=1, 2
  404.           Do k=1, estep
  405.            Do q=1, estep
  406.             l = k - q
  407.             If (l >= 1) then
  408.             F_temp(b,d,m,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*conjg(u(m,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*p)
  409.             endif
  410.            end do
  411.           end do
  412.          end do
  413.         end do
  414.        end do
  415.      End subroutine F_deriv
  416.  
  417.  
  418.  
  419.  
  420.     Subroutine G_deriv(t, p, n, u, G_temp)
  421.      implicit none
  422.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: G_temp
  423.       real(8), intent(in) :: t
  424.       complex(8), intent(in):: p
  425.       complex(8), dimension(estep), intent(in) :: n
  426.       complex(8), intent(in), dimension(2,2,estep) :: u
  427.       integer :: b,d,m,k,l,q
  428.  
  429.        Do b=1, 2
  430.         Do d=1, 2
  431.          Do m=1, 2
  432.           Do k=1, estep
  433.            Do q=1, estep
  434.             l = k + q
  435.             If (l >= 1 .and. l<=estep) then
  436.             G_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*conjg(p)
  437.             endif
  438.            end do
  439.           end do
  440.          end do
  441.         end do
  442.        end do
  443.      End subroutine G_deriv
  444.  
  445.  
  446.  
  447.  
  448.     Subroutine H_deriv(t, p, n, u, H_temp)
  449.      implicit none
  450.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: H_temp
  451.       real(8), intent(in) :: t
  452.       complex(8), intent(in):: p
  453.       complex(8), dimension(estep), intent(in) :: n
  454.       complex(8), intent(in), dimension(2,2,estep) :: u
  455.       integer :: b,d,m,k,l,q
  456.  
  457.        Do b=1, 2
  458.         Do d=1, 2
  459.          Do m=1, 2
  460.           Do k=1, estep
  461.            Do q=1, estep
  462.             l = k + q
  463.             If (l >= 1 .and. l<=estep) then
  464.             H_temp(b,d,m,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*conjg(p))
  465.             endif
  466.            end do
  467.           end do
  468.          end do
  469.         end do
  470.        end do
  471.      End subroutine H_deriv
  472.  
  473.  
  474.  
  475.     Subroutine I_deriv(t, p, n, u, I_temp)
  476.      implicit none
  477.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: I_temp
  478.       real(8), intent(in) :: t
  479.       complex(8), intent(in):: p
  480.       complex(8), dimension(estep), intent(in) :: n
  481.       complex(8), intent(in), dimension(2,2,estep) :: u
  482.       integer :: b,d,m,k,l,q
  483.  
  484.        Do b=1, 2
  485.         Do d=1, 2
  486.          Do m=1, 2
  487.           Do k=1, estep
  488.            Do q=1, estep
  489.             l = k - q
  490.             If (l >= 1) then
  491.             I_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*conjg(u(m,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*conjg(p)
  492.             endif
  493.            end do
  494.           end do
  495.          end do
  496.         end do
  497.        end do
  498.      End subroutine I_deriv
  499.  
  500.  
  501.  
  502.     Subroutine J_deriv(t, p, n, u, J_temp)
  503.      implicit none
  504.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: J_temp
  505.       real(8), intent(in) :: t
  506.       complex(8), intent(in):: p
  507.       complex(8), dimension(estep), intent(in) :: n
  508.       complex(8), intent(in), dimension(2,2,estep) :: u
  509.       integer :: b,d,m,k,l,q
  510.  
  511.        Do b=1, 2
  512.         Do d=1, 2
  513.          Do m=1, 2
  514.           Do k=1, estep
  515.            Do q=1, estep
  516.             l = k - q
  517.             If (l >= 1) then
  518.             J_temp(b,d,m,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*conjg(u(m,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*conjg(p))
  519.             endif
  520.            end do
  521.           end do
  522.          end do
  523.         end do
  524.        end do
  525.      End subroutine J_deriv
  526.  
  527.  
  528.  
  529.     Subroutine K_deriv(t, p, n, u, K_temp)
  530.      implicit none
  531.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: K_temp
  532.       real(8), intent(in) :: t
  533.       complex(8), intent(in):: p
  534.       complex(8), dimension(estep), intent(in) :: n
  535.       complex(8), intent(in), dimension(2,2,estep) :: u
  536.       integer :: b,d,m,k,l,q
  537.  
  538.        Do b=1, 2
  539.         Do d=1, 2
  540.          Do m=1, 2
  541.           Do k=1, estep
  542.            Do q=1, estep
  543.             l = k + q
  544.             If (l >= 1 .and. l<=estep) then
  545.             K_temp(b,d,m,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*p
  546.             endif
  547.            end do
  548.           end do
  549.          end do
  550.         end do
  551.        end do
  552.      End subroutine K_deriv
  553.  
  554.  
  555.  
  556.     Subroutine L_deriv(t, p, n, u, L_temp)
  557.      implicit none
  558.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: L_temp
  559.       real(8), intent(in) :: t
  560.       complex(8), intent(in):: p
  561.       complex(8), dimension(estep), intent(in) :: n
  562.       complex(8), intent(in), dimension(2,2,estep) :: u
  563.       integer :: b,d,m,k,l,q
  564.  
  565.        Do b=1, 2
  566.         Do d=1, 2
  567.          Do m=1, 2
  568.           Do k=1, estep
  569.            Do q=1, estep
  570.             l = k + q
  571.             If (l >= 1 .and. l<=estep) then
  572.             L_temp(b,d,m,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*u(m,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*p)
  573.             endif
  574.            end do
  575.           end do
  576.          end do
  577.         end do
  578.        end do
  579.      End subroutine L_deriv
  580.  
  581.  
  582. !trong subroutine n_deriv nay, ta su dung function "integral()" de tinh tich phan.
  583. ! Trong khai bao duoi day, su dung E_in, F_in...de tranh nham lan voi cac chi so dem i,k,m...
  584.  
  585.      Subroutine n_deriv(t, p, n, u, E_in, F_in, G_in, H_in, I_in, J_in, K_in, L_in, kn)
  586.       implicit none
  587.       complex(8), intent(in), dimension(2,2,2,estep,estep,estep) :: E_in, F_in, G_in, H_in, I_in, J_in, K_in, L_in
  588.       real(8), intent(in) :: t
  589.       complex(8), intent(in):: p
  590.       complex(8), dimension(estep), intent(in) :: n
  591.       complex(8), intent(in), dimension(2,2,estep) :: u
  592.       complex(8), intent(out), dimension(estep):: kn
  593.       integer :: a,b,c,d,i,m,k,l,q
  594.       complex(8), dimension(estep) :: summ
  595.  
  596.       kn=(0.0,0.0)
  597.       summ = (0.0,0.0)
  598.  
  599.         do k=1, estep
  600.          do q=1, estep
  601.           do l=1, estep
  602.            do a=1,2
  603.             do b=1,2
  604.              do c=1,2
  605.               do d=1,2
  606.                do i=1,2
  607.                 do m=1,2
  608.                  
  609.                  if (l >= abs(k-q) .and.  l <= abs(k+q) .and. (4*k*k*q*q-(k*k+q*q-l*l)**2) /= 0) then
  610.                   kn(k) = kn(k) + summ(k)
  611.                   summ(k) = -(f0*f0/(pi*pi))*                                                                    &
  612.                             ( (sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  613.                                 conjg(u(1,a,k))*u(1,c,l)*u(1,i,q)*E_in(b,d,m,k,l,q)                              &
  614.                              +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  615.                                 u(1,a,k)*conjg(u(1,c,l))*conjg(u(1,i,q))*F_in(b,d,m,k,l,q)                       &
  616.                              -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  617.                                 conjg(u(1,a,k))*u(1,c,l)*conjg(u(1,i,q))*G_in(b,d,m,k,l,q)                       &
  618.                              -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  619.                                 u(1,a,k)*conjg(u(1,c,l))*u(1,i,q)*H_in(b,d,m,k,l,q)                              &
  620.                              +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  621.                                 conjg(u(2,a,k))*u(2,c,l)*u(2,i,q)*I_in(b,d,m,k,l,q)                              &
  622.                              +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  623.                                 u(2,a,k)*conjg(u(2,c,l))*conjg(u(2,i,q))*J_in(b,d,m,k,l,q)                       &
  624.                              -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  625.                                 conjg(u(2,a,k))*u(2,c,l)*conjg(u(2,i,q))*K_in(b,d,m,k,l,q)                       &
  626.                              -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  627.                                 u(2,a,k)*conjg(u(2,c,l))*u(2,i,q)*L_in(b,d,m,k,l,q)  )
  628.                   end if
  629.                 end do
  630.                end do
  631.               end do
  632.              end do
  633.             end do
  634.            end do
  635.           end do
  636.          end do
  637.         end do
  638.      end subroutine n_deriv
  639.  
  640.  
  641.  
  642.  
  643.      subroutine p_deriv(t, p, n, u, A_in, B_in, C_in, D_in, kp)
  644.  
  645.       implicit none
  646.       complex(8), intent(in), dimension(2,2,2,estep,estep,estep) :: A_in, B_in, C_in, D_in
  647.       real(8), intent(in) :: t
  648.       complex(8), intent(in):: p
  649.       complex(8), dimension(estep), intent(in) :: n
  650.       complex(8), intent(in), dimension(2,2,estep) :: u
  651.       complex(8), intent(out):: kp
  652.       integer :: a,b,c,d,i,m,k,l,q
  653.       complex(8):: sump
  654.  
  655.       kp=(0.0,0.0)
  656.       sump = (0.0,0.0)
  657.  
  658.         do k=1, estep
  659.          do q=1, estep
  660.           do l=1, estep
  661.            do a=1,2
  662.             do b=1,2
  663.              do c=1,2
  664.               do d=1,2
  665.                do i=1,2
  666.                 do m=1,2
  667.                  
  668.                  if (l >= abs(k-q) .and.  l <= abs(k+q).and.(4*k*k*q*q-(k*k+q*q-l*l)**2) /= 0) then
  669.                   kp = kp + sump
  670.                   sump = (f0*f0/(2.0*pi**3))*                                                                &
  671.                        ((sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(k*q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  672.                            u(1,a,k)*conjg(u(1,c,l))*conjg(u(1,i,q))*A_in(b,d,m,k,l,q)                        &
  673.                        -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(k*q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  674.                            conjg(u(1,a,k))*u(1,c,l)*conjg(u(1,i,q))*B_in(b,d,m,k,l,q)                        &
  675.                        +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(k*q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  676.                            u(2,a,k)*conjg(u(2,c,l))*conjg(u(2,i,q))*C_in(b,d,m,k,l,q)                        &
  677.                        -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(k*q*l/(sqrt(real(4*k*k*q*q-(k*k+q*q-l*l)**2))))* &
  678.                            conjg(u(2,a,k))*u(2,c,l)*conjg(u(2,i,q))*D_in(b,d,m,k,l,q) )
  679.  
  680.  
  681.                   end if
  682.                 end do
  683.                end do
  684.               end do
  685.              end do
  686.             end do
  687.            end do
  688.           end do
  689.          end do
  690.         end do
  691.     kp= conjg(kp)
  692.    
  693.     end subroutine p_deriv
  694.   End module runge
  695.  
  696.  
  697. !===================================================================================================
  698.  
  699. !Bat dau chuong trinh chinh.
  700. !Chuong trinh nay se doc cac tham so dau vao: e(k) tu file. Doc cac gia tri t, h, f0 tu nguoi dung
  701. ! Sau do se xuat ra file các gia tri cua p(t), n(k,t)
  702.  
  703.   Program process
  704.     use global
  705.     use runge
  706.     implicit none
  707.  
  708.     complex(8) :: p                                
  709.     complex(8):: p_out
  710.  
  711.     complex(8) , dimension(estep) :: n            
  712.     complex(8), dimension(estep) :: n_out
  713.  
  714.     real(8) :: t  , t_max
  715.     real(8) :: t_out                        
  716.    
  717.     complex(8),  dimension(2,2,estep):: u
  718.     complex(8),  dimension(2,2,estep):: u_out  
  719.  
  720.     complex(8), dimension(2,2,2,estep,estep,estep) :: A_in, B_in, C_in, D_in, E_in, F_in,&
  721.                                                                   G_in, H_in, I_in, J_in, K_in, L_in
  722.     complex(8), dimension(2,2,2,estep, estep, estep) ::A_out, B_out, C_out, D_out, E_out, F_out, G_out, H_out,&
  723.                                                  I_out, J_out, K_out, L_out
  724.     integer :: i
  725.  
  726.  sigma_z(1,1) = (1.0,0.0)
  727.  sigma_z(1,2) = (0.0,0.0)
  728.  sigma_z(2,1) = (0.0,0.0)
  729.  sigma_z(2,2) =(-1.0,0.0)
  730.  sigma_y(1,1) = (0.0,0.0)
  731.  sigma_y(1,2) =-iu
  732.  sigma_y(2,1) = iu
  733.  sigma_y(2,2) = (0.0,0.0)
  734.  
  735.  
  736. Print *,"Nhap t0"
  737. read *, t
  738.  
  739. Print *, "Nhap h, voi h la buoc cua t"
  740. read * , h
  741.  
  742. Print *, "Nhap t_max"
  743. read *, t_max
  744.  
  745. Print *, "Nhap e(k)"
  746.   do i=1, estep
  747.     read *, e(i)
  748.   end do
  749.  
  750. p= (1.0,0.0)
  751. n=(1.0,0.0)
  752. u=(2.0,1.0)
  753. A_in=(1.0,0.0)
  754. B_in=(1.0,0.0)
  755. C_in=(1.0,0.0)
  756. D_in=(1.0,0.0)
  757. E_in=(1.0,0.0)
  758. F_in=(1.0,0.0)
  759. G_in=(1.0,0.0)
  760. H_in=(1.0,0.0)
  761. I_in=(1.0,0.0)
  762. J_in=(1.0,0.0)
  763. K_in=(1.0,0.0)
  764. L_in=(1.0,0.0)
  765.  
  766.  
  767. f0=1.0
  768.  
  769. open (unit=1, file="c:\\result1.txt")
  770.  
  771. do while(t <= t_max)
  772.  
  773.  call rk4(t, h, e, f0 ,p, n, u,A_in, B_in, C_in, D_in, E_in, F_in, G_in, H_in,I_in, J_in, K_in, L_in,&
  774.                               A_out, B_out, C_out, D_out, E_out, F_out,G_out, H_out, I_out, J_out, K_out, L_out, &
  775.                                t_out, p_out, n_out, u_out)
  776.   print *, t_out
  777.  
  778.   write (1,*)  p_out
  779.   t = t_out
  780.   p = p_out
  781.   n = n_out
  782.   u = u_out
  783.   A_in = A_out
  784.   B_in = B_out
  785.   C_in = C_out
  786.   D_in = D_out
  787.   E_in = E_out
  788.   F_in = F_out
  789.   G_in = G_out
  790.   H_in = H_out
  791.   I_in = I_out
  792.   J_in = J_out
  793.   K_in = K_out
  794.   L_in = L_out
  795.  
  796. end do
  797.  
  798. end program process
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!