API tools faq
paste
Advertisement
Guest User

Untitled

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