API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jun 15th, 2017
88
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 17.83 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
  3. ! Cac gia tri ban dau: n(t=0)=0 ; p(t=0) = 0; u(t=0) =0; t=0;
  4.  
  5.  
  6. Module global      ! Khai bao cac hang so se su dung trong chuong trinh
  7.  implicit none
  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 :: estep  ! Luoi' k (cac gia tri cua k)
  13.  complex, parameter, 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) :: ek     ! Pho gia tri nang luong theo k.
  24.  real(8) ::   f0
  25.  
  26. End module global
  27.  
  28. !============================================================================================================
  29.  
  30.  
  31. Module Runge                                                ! Giai thuat Runge-kutta gan dung cap 4.
  32.  use global
  33.  contains                                                   ! Lenh nay dung de bao lai cac subroutine con cua subroutine Runge
  34.    subroutine rk4(t, h, p, n, u,A_in, B_in, C_in, D_in, E_in, F_in, H_in,I_in, J_in, K_in, L_in,&
  35.                               A_out, B_out, C_out, D_out, E_out, F_out, H_out, I_out, J_out, K_out, L_out, &
  36.                                t_out, p_out, n_out, u_out)
  37.                      
  38.  
  39.     implicit none
  40.  
  41.     complex(8), intent(in) :: p                                ! Vi p là ham chi phu thuoc t, ko phu thuoc k.
  42.     complex(8), intent(out):: p_out
  43.  
  44.     complex(8), intent(in) , dimension(estep) :: n             ! mang n gom estep phan tu, moi phan tu ung voi mot gia tri cua k.
  45.     complex(8), intent(out), dimension(estep) :: n_out
  46.  
  47.     real(8), intent(in) :: t ,  h  
  48.     real(8), intent(out) :: t_out                        ! t la thoi gian, h la buoc nhay cua t
  49.    
  50.     complex(8), intent (in), dimension(2,2,estep):: u
  51.     complex(8), intent (out), dimension(2,2,estep):: u_out  !Ung voi moi gia tri cua k, ham U la mot matran 2x2
  52.    
  53.     complex(8), dimension(estep):: p_deriv, n_deriv, u_deriv   ! khai bao 3 phuong trinh.cua p, n ,u
  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.  
  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 n_deriv(t,p,n,u, A_in, B_in, C_in, D_in,kn1)
  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 p_deriv(t,p,n,u, E_1, F_1, G_1, H_1, I_1, J_1, K_1, L_1, kp1)
  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 n_deriv(t1,p1,n1,u1, A_1, B_1, C_1, D_1,kn2)
  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 p_deriv(t1,p1,n1,u1, E_2, F_2, G_2, H_2, I_2, J_2, K_2, L_2, kp2)
  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 n_deriv(t1, p2, n2, u2, A_2, B_2, C_2, D_2, kn3)
  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 p_deriv(t1, p2, n2, u2, E_2, F_2, G_2, H_2, I_2, J_2, K_2, L_2, kp3)
  178.  
  179.      t3 = t + h
  180.      u3 = u + h*ku3
  181.      n3 = n + h*ku3
  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 n_deriv(t3, p3, n3, u3, A_3, B_3, C_3, D_3, kn4)
  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 p_deriv(t3, p3, n3, u3, E_3, F_3, G_3, H_3, I_3, J_3, K_3, L_3, kp4)
  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.      Do i =1 , estep                                                          ! Thiet lap H, moi phan tu cua H theo k la mot matran 2x2
  254.       H(2,2,i) = (e(i) + f0*n(i))* sigma_z + iu*f0*p*p*sigma_y
  255.      end do
  256.  
  257.      Do i = 1, estep
  258.       ku(2,2,i) = -iu* matmul(H(2,2,i),u(2,2,i))                              !Day la pt A.11 , Lenh matmul() la nhan ma tran
  259.      end do
  260.  
  261.     end subroutine u_deriv
  262.  
  263.  
  264.     Subroutine A_deriv(t, p, n, u, A_temp)
  265.      implicit none
  266.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: A_temp
  267.       real(8), intent(in) :: t ,
  268.       complex(8), intent(in):: p
  269.       complex(8), dimension(estep), intent(in) :: n
  270.       complex(8), intent(in), dimension(2,2,estep) :: u
  271.       integer :: b,d,n,k,l,q
  272.  
  273.        Do b=1, 2
  274.         Do d=1, 2
  275.          Do n=1, 2
  276.           Do k=1, estep
  277.            Do q=1, estep
  278.             l = k - q
  279.             If (l >= 1) then
  280.             A_temp(b,d,n,k,l,q) = conjg(u(b,1,k))*u(d,1,l)*u(n,1,q)*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*conjg(p)
  281.             endif
  282.            end do
  283.           end do
  284.          end do
  285.         end do
  286.        end do
  287.  
  288.  
  289.  
  290.      Subroutine B_deriv(t, p, n, u, B_temp)
  291.       implicit none
  292.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: B_temp
  293.       real(8), intent(in) :: t ,
  294.       complex(8), intent(in):: p
  295.       complex(8), dimension(estep), intent(in) :: n
  296.       complex(8), intent(in), dimension(2,2,estep) :: u
  297.       integer :: b,d,n,k,l,q
  298.        
  299.        Do b=1, 2
  300.         Do d=1, 2
  301.          Do n=1, 2
  302.           Do k=1, estep
  303.            Do q=1, estep
  304.             l = k + q
  305.             If (l >= 1) then
  306.             B_temp(b,d,n,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(n,1nq)*(n(l) + n(l)*n(q)+n(l)*n(k) - n(k)*n(q))*conjg(p)
  307.             endif
  308.            end do
  309.           end do
  310.          end do
  311.         end do
  312.        end do
  313.  
  314.      End subroutine B_deriv
  315.  
  316.  
  317.     Subroutine C_deriv(t, p, n, u, C_temp)
  318.      implicit none
  319.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: C_temp
  320.       real(8), intent(in) :: t ,
  321.       complex(8), intent(in):: p
  322.       complex(8), dimension(estep), intent(in) :: n
  323.       complex(8), intent(in), dimension(2,2,estep) :: u
  324.       integer :: b,d,n,k,l,q
  325.  
  326.        Do b=1, 2
  327.         Do d=1, 2
  328.          Do n=1, 2
  329.           Do k=1, estep
  330.            Do q=1, estep
  331.             l = k - q
  332.             If (l >= 1) then
  333.             C_temp(b,d,n,k,l,q) = conjg(u(b,1,k))*u(d,1,l)*u(n,1,q)*(n(k)+n(k)*n(l)+n(k)*n(q)-n(l)*n(q))*p
  334.             endif
  335.            end do
  336.           end do
  337.          end do
  338.         end do
  339.        end do
  340.      End subroutine C_deriv
  341.  
  342.  
  343.  
  344.     Subroutine D_deriv(t, p, n, u, D_temp)
  345.      implicit none
  346.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: D_temp
  347.       real(8), intent(in) :: t ,
  348.       complex(8), intent(in):: p
  349.       complex(8), dimension(estep), intent(in) :: n
  350.       complex(8), intent(in), dimension(2,2,estep) :: u
  351.       integer :: b,d,n,k,l,q
  352.  
  353.        Do b=1, 2
  354.         Do d=1, 2
  355.          Do n=1, 2
  356.           Do k=1, estep
  357.            Do q=1, estep
  358.             l = k + q
  359.             If (l >= 1) then
  360.             D_temp(b,d,n,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(n,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*p
  361.             endif
  362.            end do
  363.           end do
  364.          end do
  365.         end do
  366.        end do
  367.      End subroutine D_deriv
  368.  
  369.  
  370.  
  371.     Subroutine E_deriv(t, p, n, u, E_temp)
  372.      implicit none
  373.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: E_temp
  374.       real(8), intent(in) :: t ,
  375.       complex(8), intent(in):: p
  376.       complex(8), dimension(estep), intent(in) :: n
  377.       complex(8), intent(in), dimension(2,2,estep) :: u
  378.       integer :: b,d,n,k,l,q
  379.  
  380.        Do b=1, 2
  381.         Do d=1, 2
  382.          Do n=1, 2
  383.           Do k=1, estep
  384.            Do q=1, estep
  385.             l = k - q
  386.             If (l >= 1) then
  387.             E_temp(b,d,n,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*conjg(u(n,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*p
  388.             endif
  389.            end do
  390.           end do
  391.          end do
  392.         end do
  393.        end do
  394.      End subroutine E_deriv
  395.  
  396.  
  397. !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.
  398.     Subroutine F_deriv(t, p, n, u, F_temp)
  399.      implicit none
  400.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: F_temp
  401.       real(8), intent(in) :: t ,
  402.       complex(8), intent(in):: p
  403.       complex(8), dimension(estep), intent(in) :: n
  404.       complex(8), intent(in), dimension(2,2,estep) :: u
  405.       integer :: b,d,n,k,l,q
  406.  
  407.        Do b=1, 2
  408.         Do d=1, 2
  409.          Do n=1, 2
  410.           Do k=1, estep
  411.            Do q=1, estep
  412.             l = k - q
  413.             If (l >= 1) then
  414.             F_temp(b,d,n,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*conjg(u(n,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*p)
  415.             endif
  416.            end do
  417.           end do
  418.          end do
  419.         end do
  420.        end do
  421.      End subroutine F_deriv
  422.  
  423.  
  424.  
  425.  
  426.     Subroutine G_deriv(t, p, n, u, G_temp)
  427.      implicit none
  428.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: G_temp
  429.       real(8), intent(in) :: t ,
  430.       complex(8), intent(in):: p
  431.       complex(8), dimension(estep), intent(in) :: n
  432.       complex(8), intent(in), dimension(2,2,estep) :: u
  433.       integer :: b,d,n,k,l,q
  434.  
  435.        Do b=1, 2
  436.         Do d=1, 2
  437.          Do n=1, 2
  438.           Do k=1, estep
  439.            Do q=1, estep
  440.             l = k + q
  441.             If (l >= 1) then
  442.             G_temp(b,d,n,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(n,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*conjg(p)
  443.             endif
  444.            end do
  445.           end do
  446.          end do
  447.         end do
  448.        end do
  449.      End subroutine G_deriv
  450.  
  451.  
  452.  
  453.  
  454.     Subroutine H_deriv(t, p, n, u, H_temp)
  455.      implicit none
  456.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: H_temp
  457.       real(8), intent(in) :: t ,
  458.       complex(8), intent(in):: p
  459.       complex(8), dimension(estep), intent(in) :: n
  460.       complex(8), intent(in), dimension(2,2,estep) :: u
  461.       integer :: b,d,n,k,l,q
  462.  
  463.        Do b=1, 2
  464.         Do d=1, 2
  465.          Do n=1, 2
  466.           Do k=1, estep
  467.            Do q=1, estep
  468.             l = k + q
  469.             If (l >= 1) then
  470.             H_temp(b,d,n,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*u(n,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*conjg(p))
  471.             endif
  472.            end do
  473.           end do
  474.          end do
  475.         end do
  476.        end do
  477.      End subroutine H_deriv
  478.  
  479.  
  480.  
  481.     Subroutine I_deriv(t, p, n, u, I_temp)
  482.      implicit none
  483.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: I_temp
  484.       real(8), intent(in) :: t ,
  485.       complex(8), intent(in):: p
  486.       complex(8), dimension(estep), intent(in) :: n
  487.       complex(8), intent(in), dimension(2,2,estep) :: u
  488.       integer :: b,d,n,k,l,q
  489.  
  490.        Do b=1, 2
  491.         Do d=1, 2
  492.          Do n=1, 2
  493.           Do k=1, estep
  494.            Do q=1, estep
  495.             l = k - q
  496.             If (l >= 1) then
  497.             I_temp(b,d,n,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*conjg(u(n,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*conjg(p)
  498.             endif
  499.            end do
  500.           end do
  501.          end do
  502.         end do
  503.        end do
  504.      End subroutine I_deriv
  505.  
  506.  
  507.  
  508.     Subroutine J_deriv(t, p, n, u, J_temp)
  509.      implicit none
  510.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: J_temp
  511.       real(8), intent(in) :: t ,
  512.       complex(8), intent(in):: p
  513.       complex(8), dimension(estep), intent(in) :: n
  514.       complex(8), intent(in), dimension(2,2,estep) :: u
  515.       integer :: b,d,n,k,l,q
  516.  
  517.        Do b=1, 2
  518.         Do d=1, 2
  519.          Do n=1, 2
  520.           Do k=1, estep
  521.            Do q=1, estep
  522.             l = k - q
  523.             If (l >= 1) then
  524.             J_temp(b,d,n,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*conjg(u(n,1,q))*(n(k)+n(k)*n(q)+n(k)*n(l)-n(l)*n(q))*conjg(p))
  525.             endif
  526.            end do
  527.           end do
  528.          end do
  529.         end do
  530.        end do
  531.      End subroutine J_deriv
  532.  
  533.  
  534.  
  535.     Subroutine K_deriv(t, p, n, u, K_temp)
  536.      implicit none
  537.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: K_temp
  538.       real(8), intent(in) :: t ,
  539.       complex(8), intent(in):: p
  540.       complex(8), dimension(estep), intent(in) :: n
  541.       complex(8), intent(in), dimension(2,2,estep) :: u
  542.       integer :: b,d,n,k,l,q
  543.  
  544.        Do b=1, 2
  545.         Do d=1, 2
  546.          Do n=1, 2
  547.           Do k=1, estep
  548.            Do q=1, estep
  549.             l = k + q
  550.             If (l >= 1) then
  551.             K_temp(b,d,n,k,l,q) = u(b,1,k)*conjg(u(d,1,l))*u(n,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*p
  552.             endif
  553.            end do
  554.           end do
  555.          end do
  556.         end do
  557.        end do
  558.      End subroutine K_deriv
  559.  
  560.  
  561.  
  562.     Subroutine L_deriv(t, p, n, u, L_temp)
  563.      implicit none
  564.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: L_temp
  565.       real(8), intent(in) :: t ,
  566.       complex(8), intent(in):: p
  567.       complex(8), dimension(estep), intent(in) :: n
  568.       complex(8), intent(in), dimension(2,2,estep) :: u
  569.       integer :: b,d,n,k,l,q
  570.  
  571.        Do b=1, 2
  572.         Do d=1, 2
  573.          Do n=1, 2
  574.           Do k=1, estep
  575.            Do q=1, estep
  576.             l = k + q
  577.             If (l >= 1) then
  578.             L_temp(b,d,n,k,l,q) = conjg(u(b,1,k)*conjg(u(d,1,l))*u(n,1,q)*(n(l)+n(l)*n(q)+n(l)*n(k)-n(k)*n(q))*p)
  579.             endif
  580.            end do
  581.           end do
  582.          end do
  583.         end do
  584.        end do
  585.      End subroutine L_deriv
  586.  
  587.  
  588. !trong subroutine n_deriv nay, ta su dung function "integral()" de tinh tich phan.
  589.      Subroutine n_deriv(t, p, n, u, A, B, C, D, kn)
  590.       implicit none
  591.       complex(8), intent(in), dimension(2,2,2,estep,estep,estep) :: A_in,B_in,C_in,D_in
  592.       real(8), intent(in) :: t ,
  593.       complex(8), intent(in):: p
  594.       complex(8), dimension(estep), intent(in) :: n
  595.       complex(8), intent(in), dimension(2,2,estep) :: u
  596.       complex(8), intent(out), dimension(step):: kn
  597.       integer :: b,d,n,k,l,q
  598.  
  599.  
  600.  
  601.  
  602.  
  603.  
  604.  
  605.  
  606.  
  607.  
  608.  
  609.  
  610.  
  611.  
  612.  
  613.  
  614.  
  615.     Subroutine p_deriv(t,p,n,u,kn)
  616.      implicit none
  617.  
  618.      real(8), intent(in) :: t
  619.      real(8), intent(in) ::p
  620.      real(8), intent(in), dimension(estep) :: n
  621.      complex(8), intent(in), dimension(2,2,estep) :: u
  622.      real(8), intent(out), dimension(estep) :: kn
  623.  
  624.      integer :: i
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!