API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jun 19th, 2017
76
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 28.75 KB | None | 0 0
raw download clone embed print report
  1.  
  2. ! Chuong trinh se tinh n(t), p(t) và u(t)
  3.  
  4. ! tham so dau vao: e(k), f0   (cac tham so nay se dc nhap vao, hoac viet mot truong chinh tinh toan rieng.
  5.  
  6. ! Cac gia tri ban dau: n(t=0)=0 ; p(t=0) = 10^-6; u(t=0) =1.; t0=0;
  7.  
  8.  
  9.  
  10. Module global      ! Khai bao cac hang so se su dung trong chuong trinh
  11.  
  12.  implicit none
  13.  
  14.  save
  15.  
  16.  real, parameter :: pi= 3.141593
  17.  
  18.  real, parameter :: hbar = 658.d-3
  19.  
  20.  complex, parameter :: iu = (0.,1.) , ui = (1.,0.)
  21.  
  22.  real(8) :: t0
  23.  
  24.  integer, parameter, public:: estep = 10  ! Luoi' k (cac gia tri cua k)
  25.  
  26.  complex, dimension(2,2) :: sigma_z, sigma_y  ! khai bao 2 ma tran sigma
  27.  
  28.  
  29.  
  30. ! Pho gia tri nang luong theo k.
  31.  
  32.  real(8) ::    h
  33.  
  34. real(8),  dimension(estep) ::E
  35.  
  36.  real, parameter :: Lc = 10000.0                                           ! don vi nm
  37.  
  38.  real, parameter :: E0 = 1.616d3 , Ec0 = 0.1 + 1.616d3                                       ! don vi meV
  39.  
  40.  real, parameter :: hbar_omega = 4.3                                            ! don vi meV
  41.  
  42.  real(8), parameter :: m_e = 0.067* 5.7d-3 , m_h = 0.45 *5.7d-3
  43.  
  44.  real(8), parameter :: tau_c = 2.0    , tau_x= 50.0                                                 ! don vi ps
  45.  
  46.  real(8), parameter :: S = 10.0**8                                                 !nm^2
  47.  
  48. End module global
  49.  
  50.  
  51.  
  52. !============================================================================================================
  53.  
  54.  
  55.  
  56.  
  57.  
  58. Module Runge                                                ! Runge-kutta gan dung bac 4.
  59.  
  60.  use global
  61.  
  62.  implicit none
  63.  
  64.  contains                                                   ! Lenh nay dung de bao lai cac subroutine con cua subroutine Runge
  65.  
  66.    subroutine rk4(t, h, e, 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,&
  67.  
  68.                               A_out, B_out, C_out, D_out, E_out, F_out, G_out, H_out, I_out, J_out, K_out, L_out, &
  69.  
  70.                                t_out, p_out, n_out, u_out)
  71.  
  72.                      
  73.  
  74.  
  75.  
  76.     implicit none
  77.  
  78.     real(8), intent(in), dimension(estep) ::e
  79.  
  80.     complex(8), intent(in) :: p                                ! Vi p là ham chi phu thuoc t, ko phu thuoc k.
  81.  
  82.     complex(8), intent(out):: p_out
  83.  
  84.  
  85.  
  86.     complex(8), intent(in) , dimension(estep) :: n             ! mang n gom estep phan tu, moi phan tu ung voi mot gia tri cua k.
  87.  
  88.     complex(8), intent(out), dimension(estep) :: n_out
  89.  
  90.  
  91.  
  92.     real(8), intent(in) :: t ,  h  
  93.  
  94.     real(8), intent(out) :: t_out                        ! t la thoi gian, h la buoc nhay cua t
  95.  
  96.    
  97.  
  98.     complex(8), intent (in), dimension(2,2,estep):: u
  99.  
  100.     complex(8), intent (out), dimension(2,2,estep):: u_out  !Ung voi moi gia tri cua k, ham U la mot matran 2x2
  101.  
  102.  
  103.  
  104.     complex(8), intent(in),dimension(2,2,2,estep,estep,estep) :: A_in, B_in, C_in, D_in, E_in, F_in,&
  105.  
  106.                                                                   G_in, H_in, I_in, J_in, K_in, L_in
  107.  
  108.     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, &
  109.  
  110.                                                  I_out, J_out, K_out, L_out
  111.  
  112.  
  113.  
  114.  
  115.  
  116.     ! Khai bao cac he so tam thoi trong subroutine:
  117.  
  118.     real(8) :: t1, t3
  119.  
  120.     complex(8) :: kp1, kp2, kp3, kp4, p1, p2, p3
  121.  
  122.     complex(8), dimension(estep) :: kn1, kn2, kn3, kn4, n1, n2, n3
  123.  
  124.     complex(8), dimension(2,2,estep):: ku1, ku2, ku3, ku4, u1, u2, u3
  125.  
  126.     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,  &
  127.  
  128.                                                          C_temp1, C_temp2, C_temp3, C_temp4, D_temp1, D_temp2, D_temp3, D_temp4  
  129.  
  130.  
  131.  
  132.  
  133.  
  134.     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, &
  135.  
  136.                                                    G_temp1, G_temp2, G_temp3, G_temp4, H_temp1, H_temp2, H_temp3, H_temp4, &
  137.  
  138.                                                    I_temp1, I_temp2, I_temp3, I_temp4, J_temp1, J_temp2, J_temp3, J_temp4, &
  139.  
  140.                                                    K_temp1, K_temp2, K_temp3, K_temp4, L_temp1, L_temp2, L_temp3, L_temp4
  141.  
  142.    
  143.  
  144.     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, &
  145.  
  146.                                                      E_1, E_2, E_3, F_1, F_2, F_3, G_1, G_2, G_3, H_1, H_2, H_3, &
  147.  
  148.                                                      I_1, I_2, I_3, J_1, J_2, J_3, K_1, K_2, K_3, L_1, L_2, L_3                                                      
  149.  
  150.              
  151.  
  152.  
  153.  
  154.     call u_deriv(t,e,p,n,u,ku1)
  155.  
  156.  
  157.  
  158.     Call A_deriv(t,p,n,u,A_temp1)
  159.  
  160.     call B_deriv(t,p,n,u,B_temp1)
  161.  
  162.     call C_deriv(t,p,n,u,C_temp1)
  163.  
  164.     call D_deriv(t,p,n,u,D_temp1)
  165.  
  166.  
  167.  
  168.     call p_deriv(t,p,n,u, A_in, B_in, C_in, D_in,kp1)
  169.  
  170.  
  171.  
  172.     call E_deriv(t,p,n,u,E_temp1)
  173.  
  174.     call F_deriv(t,p,n,u,F_temp1)
  175.  
  176.     call G_deriv(t,p,n,u,G_temp1)
  177.  
  178.     call H_deriv(t,p,n,u,H_temp1)
  179.  
  180.     call I_deriv(t,p,n,u,I_temp1)
  181.  
  182.     call J_deriv(t,p,n,u,J_temp1)
  183.  
  184.     call K_deriv(t,p,n,u,K_temp1)
  185.  
  186.     call L_deriv(t,p,n,u,L_temp1)
  187.  
  188.  
  189.  
  190.     call n_deriv(t,p,n,u, E_1, F_1, G_1, H_1, I_1, J_1, K_1, L_1, kn1)
  191.  
  192.  
  193.  
  194.       t1  = t+ h/2.0
  195.  
  196.       u1 = u + h*ku1/2.0
  197.  
  198.       n1 = n + h*kn1/2.0
  199.  
  200.       p1 = p + h*kp1/2.0
  201.  
  202.  
  203.  
  204.       A_1 = A_in + h*A_temp1/2.0
  205.  
  206.       B_1 = B_in + h*B_temp1/2.0
  207.  
  208.       C_1 = C_in + h*C_temp1/2.0
  209.  
  210.       D_1 = D_in + h*D_temp1/2.0
  211.  
  212.       E_1 = E_in + h*E_temp1/2.0
  213.  
  214.       F_1 = F_in + h*F_temp1/2.0
  215.  
  216.       G_1 = G_in + h*G_temp1/2.0
  217.  
  218.       H_1 = H_in + h*H_temp1/2.0
  219.  
  220.       I_1 = I_in + h*I_temp1/2.0
  221.  
  222.       J_1 = J_in + h*J_temp1/2.0
  223.  
  224.       K_1 = K_in + h*K_temp1/2.0
  225.  
  226.       L_1 = L_in + h*L_temp1/2.0
  227.  
  228.  
  229.  
  230.  
  231.  
  232.  
  233.  
  234.     call u_deriv(t1, e, p1, n1, u1, ku2)
  235.  
  236.  
  237.  
  238.     Call A_deriv(t1,p1,n1,u1,A_temp2)
  239.  
  240.     call B_deriv(t1,p1,n1,u1,B_temp2)
  241.  
  242.     call C_deriv(t1,p1,n1,u1,C_temp2)
  243.  
  244.     call D_deriv(t1,p1,n1,u1,D_temp2)
  245.  
  246.  
  247.  
  248.     call p_deriv(t1,p1,n1,u1, A_1, B_1, C_1, D_1,kp2)
  249.  
  250.  
  251.  
  252.     call E_deriv(t1,p1,n1,u1,E_temp2)
  253.  
  254.     call F_deriv(t1,p1,n1,u1,F_temp2)
  255.  
  256.     call G_deriv(t1,p1,n1,u1,G_temp2)
  257.  
  258.     call H_deriv(t1,p1,n1,u1,H_temp2)
  259.  
  260.     call I_deriv(t1,p1,n1,u1,I_temp2)
  261.  
  262.     call J_deriv(t1,p1,n1,u1,J_temp2)
  263.  
  264.     call K_deriv(t1,p1,n1,u1,K_temp2)
  265.  
  266.     call L_deriv(t1,p1,n1,u1,L_temp2)
  267.  
  268.  
  269.  
  270.     call n_deriv(t1,p1,n1,u1, E_2, F_2, G_2, H_2, I_2, J_2, K_2, L_2, kn2)
  271.  
  272.      
  273.  
  274.       u2 = u + h*ku2/2.0
  275.  
  276.       n2 = n + h*kn2/2.0
  277.  
  278.       p2 = p + h*kp2/2.0
  279.  
  280.  
  281.  
  282.       A_2 = A_in + h*A_temp2/2.0
  283.  
  284.       B_2 = B_in + h*B_temp2/2.0
  285.  
  286.       C_2 = C_in + h*C_temp2/2.0
  287.  
  288.       D_2 = D_in + h*D_temp2/2.0
  289.  
  290.       E_2 = E_in + h*E_temp2/2.0
  291.  
  292.       F_2 = F_in + h*F_temp2/2.0
  293.  
  294.       G_2 = G_in + h*G_temp2/2.0
  295.  
  296.       H_2 = H_in + h*H_temp2/2.0
  297.  
  298.       I_2 = I_in + h*I_temp2/2.0
  299.  
  300.       J_2 = J_in + h*J_temp2/2.0
  301.  
  302.       K_2 = K_in + h*K_temp2/2.0
  303.  
  304.       L_2 = L_in + h*L_temp2/2.0     
  305.  
  306.      
  307.  
  308.  
  309.  
  310.  
  311.  
  312.     call u_deriv(t1, e, p2, n2, u2, ku3)
  313.  
  314.  
  315.  
  316.     Call A_deriv(t1,p2,n2,u2,A_temp3)
  317.  
  318.     call B_deriv(t1,p2,n2,u2,B_temp3)
  319.  
  320.     call C_deriv(t1,p2,n2,u2,C_temp3)
  321.  
  322.     call D_deriv(t1,p2,n2,u2,D_temp3)
  323.  
  324.  
  325.  
  326.     call p_deriv(t1, p2, n2, u2, A_2, B_2, C_2, D_2, kp3)
  327.  
  328.  
  329.  
  330.     call E_deriv(t1,p2,n2,u2,E_temp3)
  331.  
  332.     call F_deriv(t1,p2,n2,u2,F_temp3)
  333.  
  334.     call G_deriv(t1,p2,n2,u2,G_temp3)
  335.  
  336.     call H_deriv(t1,p2,n2,u2,H_temp3)
  337.  
  338.     call I_deriv(t1,p2,n2,u2,I_temp3)
  339.  
  340.     call J_deriv(t1,p2,n2,u2,J_temp3)
  341.  
  342.     call K_deriv(t1,p2,n2,u2,K_temp3)
  343.  
  344.     call L_deriv(t1,p2,n2,u2,L_temp3)
  345.  
  346.  
  347.  
  348.     call n_deriv(t1, p2, n2, u2, E_2, F_2, G_2, H_2, I_2, J_2, K_2, L_2, kn3)
  349.  
  350.  
  351.  
  352.      t3 = t + h
  353.  
  354.      u3 = u + h*ku3
  355.  
  356.      n3 = n + h*kn3
  357.  
  358.      p3 = p + h*kp3
  359.  
  360.      A_3 = A_in + h*A_temp3
  361.  
  362.      B_3 = B_in + h*B_temp3
  363.  
  364.      C_3 = C_in + h*C_temp3
  365.  
  366.      D_3 = D_in + h*D_temp3
  367.  
  368.      E_3 = E_in + h*E_temp3
  369.  
  370.      F_3 = F_in + h*F_temp3
  371.  
  372.      G_3 = G_in + h*G_temp3
  373.  
  374.      H_3 = H_in + h*H_temp3
  375.  
  376.      I_3 = I_in + h*I_temp3
  377.  
  378.      J_3 = J_in + h*J_temp3
  379.  
  380.      K_3 = K_in + h*K_temp3
  381.  
  382.      L_3 = L_in + h*L_temp3
  383.  
  384.  
  385.  
  386.      
  387.  
  388.  
  389.  
  390.     call u_deriv(t3, e, p3, n3, u3, ku4)
  391.  
  392.  
  393.  
  394.     Call A_deriv(t3,p3,n3,u3,A_temp4)
  395.  
  396.     call B_deriv(t3,p3,n3,u3,B_temp4)
  397.  
  398.     call C_deriv(t3,p3,n3,u3,C_temp4)
  399.  
  400.     call D_deriv(t3,p3,n3,u3,D_temp4)
  401.  
  402.  
  403.  
  404.     call p_deriv(t3, p3, n3, u3, A_3, B_3, C_3, D_3, kp4)
  405.  
  406.  
  407.  
  408.     call E_deriv(t3,p3,n3,u3,E_temp4)
  409.  
  410.     call F_deriv(t3,p3,n3,u3,F_temp4)
  411.  
  412.     call G_deriv(t3,p3,n3,u3,G_temp4)
  413.  
  414.     call H_deriv(t3,p3,n3,u3,H_temp4)
  415.  
  416.     call I_deriv(t3,p3,n3,u3,I_temp4)
  417.  
  418.     call J_deriv(t3,p3,n3,u3,J_temp4)
  419.  
  420.     call K_deriv(t3,p3,n3,u3,K_temp4)
  421.  
  422.     call L_deriv(t3,p3,n3,u3,L_temp4)
  423.  
  424.  
  425.  
  426.     call n_deriv(t3, p3, n3, u3, E_3, F_3, G_3, H_3, I_3, J_3, K_3, L_3, kn4)
  427.  
  428.  
  429.  
  430.      t_out = t3
  431.  
  432.      u_out = u + h*(ku1 +ku4 + 2.0*(ku2 + ku3))/6.0
  433.  
  434.      n_out = n + h*(kn1 +kn4 + 2.0*(kn2 + kn3))/6.0
  435.  
  436.      p_out = p + h*(kp1 +kp4 + 2.0*(kp2 + kp3))/6.0
  437.  
  438.  
  439.  
  440.      A_out = A_in + h*(A_temp1 + A_temp4 +2.0*(A_temp2 + A_temp3))/6.0
  441.  
  442.      B_out = B_in + h*(B_temp1 + B_temp4 +2.0*(B_temp2 + B_temp3))/6.0
  443.  
  444.      C_out = C_in + h*(C_temp1 + C_temp4 +2.0*(C_temp2 + C_temp3))/6.0
  445.  
  446.      D_out = D_in + h*(D_temp1 + D_temp4 +2.0*(D_temp2 + D_temp3))/6.0
  447.  
  448.      E_out = E_in + h*(E_temp1 + E_temp4 +2.0*(E_temp2 + E_temp3))/6.0
  449.  
  450.      F_out = F_in + h*(F_temp1 + F_temp4 +2.0*(F_temp2 + F_temp3))/6.0
  451.  
  452.      G_out = G_in + h*(G_temp1 + G_temp4 +2.0*(G_temp2 + G_temp3))/6.0
  453.  
  454.      H_out = H_in + h*(H_temp1 + H_temp4 +2.0*(H_temp2 + H_temp3))/6.0
  455.  
  456.      I_out = I_in + h*(I_temp1 + I_temp4 +2.0*(I_temp2 + I_temp3))/6.0
  457.  
  458.      J_out = J_in + h*(J_temp1 + J_temp4 +2.0*(J_temp2 + J_temp3))/6.0
  459.  
  460.      K_out = K_in + h*(K_temp1 + K_temp4 +2.0*(K_temp2 + K_temp3))/6.0
  461.  
  462.      L_out = L_in + h*(L_temp1 + L_temp4 +2.0*(L_temp2 + L_temp3))/6.0
  463.  
  464.  
  465.  
  466.     end subroutine rk4
  467.  
  468.  
  469.  
  470.  
  471.  
  472.    Subroutine u_deriv(t,e,p,n,u,ku)
  473.  
  474.  
  475.  
  476.     implicit none
  477.  
  478.  
  479.  
  480.     real(8), intent(in) :: t
  481.  
  482.     complex(8), intent(in) ::p
  483.  
  484.     complex(8), intent(in), dimension(estep) :: n
  485.  
  486.     complex(8), intent(in), dimension(2,2,estep) :: u
  487.  
  488.     complex(8), intent(out), dimension(2,2,estep) ::ku
  489.  
  490.     real(8), intent(in), dimension(estep) ::e
  491.  
  492.     real(8), external :: f0
  493.  
  494.     integer :: i
  495.  
  496.  
  497.  
  498.     complex(8), dimension(2,2, estep) :: H                         !H la hamiltonian
  499.  
  500.     complex(8), dimension(2,2) ::H_temp , u_temp , ku_temp    
  501.  
  502.                              
  503.  
  504.      Do i =1 , estep                                             ! Thiet lap H, moi phan tu cua H theo k la mot matran 2x2
  505.  
  506.       H(:,:,i) = (e(i) + f0(0,0)*n(1))* sigma_z + iu*f0(0,0)*p*p*sigma_y               !tam thoi su dung n(1) thay cho n(0)
  507.  
  508.      end do
  509.  
  510.  
  511.  
  512.      Do i = 1, estep
  513.  
  514.  
  515.  
  516.       ku(:,:,i) = -iu* matmul(H(:,:,i),u(:,:,i))                 !Day la pt A.11 , Lenh matmul() la nhan ma tran
  517.  
  518.      end do
  519.  
  520.  
  521.  
  522.     end subroutine u_deriv
  523.  
  524.  
  525.  
  526.  
  527.  
  528.     Subroutine A_deriv(t, p, n, u, A_temp)
  529.  
  530.      implicit none
  531.  
  532.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: A_temp
  533.  
  534.       real(8), intent(in) :: t
  535.  
  536.       complex(8), intent(in):: p
  537.  
  538.       complex(8), dimension(estep), intent(in) :: n
  539.  
  540.       complex(8), intent(in), dimension(2,2,estep) :: u
  541.  
  542.       integer :: b,d,m,k,l,q
  543.  
  544.  
  545.  
  546.        Do b=1, 2
  547.  
  548.         Do d=1, 2
  549.  
  550.          Do m=1, 2
  551.  
  552.           Do k=1, estep
  553.  
  554.            Do q=1, estep
  555.  
  556.             l = k - q
  557.  
  558.             If (l >= 1) then
  559.  
  560.             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)
  561.  
  562.             endif
  563.  
  564.            end do
  565.  
  566.           end do
  567.  
  568.          end do
  569.  
  570.         end do
  571.  
  572.        end do
  573.  
  574.      End subroutine A_deriv
  575.  
  576.  
  577.  
  578.  
  579.  
  580.  
  581.  
  582.      Subroutine B_deriv(t, p, n, u, B_temp)
  583.  
  584.       implicit none
  585.  
  586.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: B_temp
  587.  
  588.       real(8), intent(in) :: t
  589.  
  590.       complex(8), intent(in):: p
  591.  
  592.       complex(8), dimension(estep), intent(in) :: n
  593.  
  594.       complex(8), intent(in), dimension(2,2,estep) :: u
  595.  
  596.       integer :: b,d,m,k,l,q
  597.  
  598.        
  599.  
  600.        Do b=1, 2
  601.  
  602.         Do d=1, 2
  603.  
  604.          Do m=1, 2
  605.  
  606.           Do k=1, estep
  607.  
  608.            Do q=1, estep
  609.  
  610.             l = k + q
  611.  
  612.             If (l >= 1 .and. l <= estep) then
  613.  
  614.             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)
  615.  
  616.             endif
  617.  
  618.            end do
  619.  
  620.           end do
  621.  
  622.          end do
  623.  
  624.         end do
  625.  
  626.        end do
  627.  
  628.      End subroutine B_deriv
  629.  
  630.  
  631.  
  632.  
  633.  
  634.     Subroutine C_deriv(t, p, n, u, C_temp)
  635.  
  636.      implicit none
  637.  
  638.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: C_temp
  639.  
  640.       real(8), intent(in) :: t
  641.  
  642.       complex(8), intent(in):: p
  643.  
  644.       complex(8), dimension(estep), intent(in) :: n
  645.  
  646.       complex(8), intent(in), dimension(2,2,estep) :: u
  647.  
  648.       integer :: b,d,m,k,l,q
  649.  
  650.  
  651.  
  652.        Do b=1, 2
  653.  
  654.         Do d=1, 2
  655.  
  656.          Do m=1, 2
  657.  
  658.           Do k=1, estep
  659.  
  660.            Do q=1, estep
  661.  
  662.             l = k - q
  663.  
  664.             If (l >= 1) then
  665.  
  666.             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
  667.  
  668.             endif
  669.  
  670.            end do
  671.  
  672.           end do
  673.  
  674.          end do
  675.  
  676.         end do
  677.  
  678.        end do
  679.  
  680.      End subroutine C_deriv
  681.  
  682.  
  683.  
  684.  
  685.  
  686.  
  687.  
  688.     Subroutine D_deriv(t, p, n, u, D_temp)
  689.  
  690.      implicit none
  691.  
  692.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: D_temp
  693.  
  694.       real(8), intent(in) :: t
  695.  
  696.       complex(8), intent(in):: p
  697.  
  698.       complex(8), dimension(estep), intent(in) :: n
  699.  
  700.       complex(8), intent(in), dimension(2,2,estep) :: u
  701.  
  702.       integer :: b,d,m,k,l,q
  703.  
  704.  
  705.  
  706.        Do b=1, 2
  707.  
  708.         Do d=1, 2
  709.  
  710.          Do m=1, 2
  711.  
  712.           Do k=1, estep
  713.  
  714.            Do q=1, estep
  715.  
  716.             l = k + q
  717.  
  718.             If (l >= 1 .and. l <= estep) then
  719.  
  720.             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
  721.  
  722.             endif
  723.  
  724.            end do
  725.  
  726.           end do
  727.  
  728.          end do
  729.  
  730.         end do
  731.  
  732.        end do
  733.  
  734.      End subroutine D_deriv
  735.  
  736.  
  737.  
  738.  
  739.  
  740.  
  741.  
  742.     Subroutine E_deriv(t, p, n, u, E_temp)
  743.  
  744.      implicit none
  745.  
  746.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: E_temp
  747.  
  748.       real(8), intent(in) :: t
  749.  
  750.       complex(8), intent(in):: p
  751.  
  752.       complex(8), dimension(estep), intent(in) :: n
  753.  
  754.       complex(8), intent(in), dimension(2,2,estep) :: u
  755.  
  756.       integer :: b,d,m,k,l,q
  757.  
  758.  
  759.  
  760.        Do b=1, 2
  761.  
  762.         Do d=1, 2
  763.  
  764.          Do m=1, 2
  765.  
  766.           Do k=1, estep
  767.  
  768.            Do q=1, estep
  769.  
  770.             l = k - q
  771.  
  772.             If (l >= 1) then
  773.  
  774.             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
  775.  
  776.             endif
  777.  
  778.            end do
  779.  
  780.           end do
  781.  
  782.          end do
  783.  
  784.         end do
  785.  
  786.        end do
  787.  
  788.      End subroutine E_deriv
  789.  
  790.  
  791.  
  792.  
  793.  
  794. !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.
  795.  
  796.     Subroutine F_deriv(t, p, n, u, F_temp)
  797.  
  798.      implicit none
  799.  
  800.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: F_temp
  801.  
  802.       real(8), intent(in) :: t
  803.  
  804.       complex(8), intent(in):: p
  805.  
  806.       complex(8), dimension(estep), intent(in) :: n
  807.  
  808.       complex(8), intent(in), dimension(2,2,estep) :: u
  809.  
  810.       integer :: b,d,m,k,l,q
  811.  
  812.  
  813.  
  814.        Do b=1, 2
  815.  
  816.         Do d=1, 2
  817.  
  818.          Do m=1, 2
  819.  
  820.           Do k=1, estep
  821.  
  822.            Do q=1, estep
  823.  
  824.             l = k - q
  825.  
  826.             If (l >= 1) then
  827.  
  828.             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)
  829.  
  830.             endif
  831.  
  832.            end do
  833.  
  834.           end do
  835.  
  836.          end do
  837.  
  838.         end do
  839.  
  840.        end do
  841.  
  842.      End subroutine F_deriv
  843.  
  844.  
  845.  
  846.  
  847.  
  848.  
  849.  
  850.  
  851.  
  852.     Subroutine G_deriv(t, p, n, u, G_temp)
  853.  
  854.      implicit none
  855.  
  856.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: G_temp
  857.  
  858.       real(8), intent(in) :: t
  859.  
  860.       complex(8), intent(in):: p
  861.  
  862.       complex(8), dimension(estep), intent(in) :: n
  863.  
  864.       complex(8), intent(in), dimension(2,2,estep) :: u
  865.  
  866.       integer :: b,d,m,k,l,q
  867.  
  868.  
  869.  
  870.        Do b=1, 2
  871.  
  872.         Do d=1, 2
  873.  
  874.          Do m=1, 2
  875.  
  876.           Do k=1, estep
  877.  
  878.            Do q=1, estep
  879.  
  880.             l = k + q
  881.  
  882.             If (l >= 1 .and. l<=estep) then
  883.  
  884.             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)
  885.  
  886.             endif
  887.  
  888.            end do
  889.  
  890.           end do
  891.  
  892.          end do
  893.  
  894.         end do
  895.  
  896.        end do
  897.  
  898.      End subroutine G_deriv
  899.  
  900.  
  901.  
  902.  
  903.  
  904.  
  905.  
  906.  
  907.  
  908.     Subroutine H_deriv(t, p, n, u, H_temp)
  909.  
  910.      implicit none
  911.  
  912.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: H_temp
  913.  
  914.       real(8), intent(in) :: t
  915.  
  916.       complex(8), intent(in):: p
  917.  
  918.       complex(8), dimension(estep), intent(in) :: n
  919.  
  920.       complex(8), intent(in), dimension(2,2,estep) :: u
  921.  
  922.       integer :: b,d,m,k,l,q
  923.  
  924.  
  925.  
  926.        Do b=1, 2
  927.  
  928.         Do d=1, 2
  929.  
  930.          Do m=1, 2
  931.  
  932.           Do k=1, estep
  933.  
  934.            Do q=1, estep
  935.  
  936.             l = k + q
  937.  
  938.             If (l >= 1 .and. l<=estep) then
  939.  
  940.             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))
  941.  
  942.             endif
  943.  
  944.            end do
  945.  
  946.           end do
  947.  
  948.          end do
  949.  
  950.         end do
  951.  
  952.        end do
  953.  
  954.      End subroutine H_deriv
  955.  
  956.  
  957.  
  958.  
  959.  
  960.  
  961.  
  962.     Subroutine I_deriv(t, p, n, u, I_temp)
  963.  
  964.      implicit none
  965.  
  966.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: I_temp
  967.  
  968.       real(8), intent(in) :: t
  969.  
  970.       complex(8), intent(in):: p
  971.  
  972.       complex(8), dimension(estep), intent(in) :: n
  973.  
  974.       complex(8), intent(in), dimension(2,2,estep) :: u
  975.  
  976.       integer :: b,d,m,k,l,q
  977.  
  978.  
  979.  
  980.        Do b=1, 2
  981.  
  982.         Do d=1, 2
  983.  
  984.          Do m=1, 2
  985.  
  986.           Do k=1, estep
  987.  
  988.            Do q=1, estep
  989.  
  990.             l = k - q
  991.  
  992.             If (l >= 1) then
  993.  
  994.             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)
  995.  
  996.             endif
  997.  
  998.            end do
  999.  
  1000.           end do
  1001.  
  1002.          end do
  1003.  
  1004.         end do
  1005.  
  1006.        end do
  1007.  
  1008.      End subroutine I_deriv
  1009.  
  1010.  
  1011.  
  1012.  
  1013.  
  1014.  
  1015.  
  1016.     Subroutine J_deriv(t, p, n, u, J_temp)
  1017.  
  1018.      implicit none
  1019.  
  1020.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: J_temp
  1021.  
  1022.       real(8), intent(in) :: t
  1023.  
  1024.       complex(8), intent(in):: p
  1025.  
  1026.       complex(8), dimension(estep), intent(in) :: n
  1027.  
  1028.       complex(8), intent(in), dimension(2,2,estep) :: u
  1029.  
  1030.       integer :: b,d,m,k,l,q
  1031.  
  1032.  
  1033.  
  1034.        Do b=1, 2
  1035.  
  1036.         Do d=1, 2
  1037.  
  1038.          Do m=1, 2
  1039.  
  1040.           Do k=1, estep
  1041.  
  1042.            Do q=1, estep
  1043.  
  1044.             l = k - q
  1045.  
  1046.             If (l >= 1) then
  1047.  
  1048.             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))
  1049.  
  1050.             endif
  1051.  
  1052.            end do
  1053.  
  1054.           end do
  1055.  
  1056.          end do
  1057.  
  1058.         end do
  1059.  
  1060.        end do
  1061.  
  1062.      End subroutine J_deriv
  1063.  
  1064.  
  1065.  
  1066.  
  1067.  
  1068.  
  1069.  
  1070.     Subroutine K_deriv(t, p, n, u, K_temp)
  1071.  
  1072.      implicit none
  1073.  
  1074.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: K_temp
  1075.  
  1076.       real(8), intent(in) :: t
  1077.  
  1078.       complex(8), intent(in):: p
  1079.  
  1080.       complex(8), dimension(estep), intent(in) :: n
  1081.  
  1082.       complex(8), intent(in), dimension(2,2,estep) :: u
  1083.  
  1084.       integer :: b,d,m,k,l,q
  1085.  
  1086.  
  1087.  
  1088.        Do b=1, 2
  1089.  
  1090.         Do d=1, 2
  1091.  
  1092.          Do m=1, 2
  1093.  
  1094.           Do k=1, estep
  1095.  
  1096.            Do q=1, estep
  1097.  
  1098.             l = k + q
  1099.  
  1100.             If (l >= 1 .and. l<=estep) then
  1101.  
  1102.             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
  1103.  
  1104.             endif
  1105.  
  1106.            end do
  1107.  
  1108.           end do
  1109.  
  1110.          end do
  1111.  
  1112.         end do
  1113.  
  1114.        end do
  1115.  
  1116.      End subroutine K_deriv
  1117.  
  1118.  
  1119.  
  1120.  
  1121.  
  1122.  
  1123.  
  1124.     Subroutine L_deriv(t, p, n, u, L_temp)
  1125.  
  1126.      implicit none
  1127.  
  1128.       complex(8), intent(out), dimension(2,2,2,estep,estep,estep) :: L_temp
  1129.  
  1130.       real(8), intent(in) :: t
  1131.  
  1132.       complex(8), intent(in):: p
  1133.  
  1134.       complex(8), dimension(estep), intent(in) :: n
  1135.  
  1136.       complex(8), intent(in), dimension(2,2,estep) :: u
  1137.  
  1138.       integer :: b,d,m,k,l,q
  1139.  
  1140.  
  1141.  
  1142.        Do b=1, 2
  1143.  
  1144.         Do d=1, 2
  1145.  
  1146.          Do m=1, 2
  1147.  
  1148.           Do k=1, estep
  1149.  
  1150.            Do q=1, estep
  1151.  
  1152.             l = k + q
  1153.  
  1154.             If (l >= 1 .and. l<=estep) then
  1155.  
  1156.             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)
  1157.  
  1158.             endif
  1159.  
  1160.            end do
  1161.  
  1162.           end do
  1163.  
  1164.          end do
  1165.  
  1166.         end do
  1167.  
  1168.        end do
  1169.  
  1170.      End subroutine L_deriv
  1171.  
  1172.  
  1173.  
  1174.  
  1175.  
  1176.  
  1177.  
  1178. ! Trong khai bao duoi day, su dung E_in, F_in...de tranh nham lan voi cac chi so dem i,k,m...
  1179.  
  1180.  
  1181.  
  1182.      Subroutine n_deriv(t, p, n, u, E_in, F_in, G_in, H_in, I_in, J_in, K_in, L_in, kn)
  1183.  
  1184.       implicit none
  1185.  
  1186.       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
  1187.  
  1188.       real(8), intent(in) :: t
  1189.  
  1190.       complex(8), intent(in):: p
  1191.  
  1192.       complex(8), dimension(estep), intent(in) :: n
  1193.  
  1194.       complex(8), intent(in), dimension(2,2,estep) :: u
  1195.  
  1196.       complex(8), intent(out), dimension(estep):: kn
  1197.  
  1198.       integer :: a,b,c,d,i,m,k,l,q
  1199.  
  1200.       complex(8), dimension(estep) :: summ
  1201.  
  1202.       real(8), external :: f0 , integral    , pump_n                                    ! khai bao ham f0
  1203.  
  1204.  
  1205.  
  1206.       kn=(0.0,0.0)
  1207.  
  1208.       summ = (0.0,0.0)
  1209.  
  1210.              
  1211.  
  1212.         do k=1, estep
  1213.  
  1214.          do q=1, estep
  1215.  
  1216.           do l=1, estep
  1217.  
  1218.            do a=1,2
  1219.  
  1220.             do b=1,2
  1221.  
  1222.              do c=1,2
  1223.  
  1224.               do d=1,2
  1225.  
  1226.                do i=1,2
  1227.  
  1228.                 do m=1,2
  1229.  
  1230.                  
  1231.  
  1232.                  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
  1233.  
  1234.                  
  1235.  
  1236.                   summ(k) = -(f0(k,q)*f0(k,q)*S/(2.0*pi*pi))*                                                                    &
  1237.  
  1238.                             ( (sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  1239.  
  1240.                                 conjg(u(1,a,k))*u(1,c,l)*u(1,i,q)*E_in(b,d,m,k,l,q)*p                              &
  1241.  
  1242.                              +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  1243.  
  1244.                                 u(1,a,k)*conjg(u(1,c,l))*conjg(u(1,i,q))*F_in(b,d,m,k,l,q)*conjg(p)                      &
  1245.  
  1246.                              -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  1247.  
  1248.                                 conjg(u(1,a,k))*u(1,c,l)*conjg(u(1,i,q))*G_in(b,d,m,k,l,q) *conjg(p)                     &
  1249.  
  1250.                              -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  1251.  
  1252.                                 u(1,a,k)*conjg(u(1,c,l))*u(1,i,q)*H_in(b,d,m,k,l,q)*p                             &
  1253.  
  1254.                              +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  1255.  
  1256.                                 conjg(u(2,a,k))*u(2,c,l)*u(2,i,q)*I_in(b,d,m,k,l,q)*conjg(p)                          &
  1257.  
  1258.                              +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  1259.  
  1260.                                 u(2,a,k)*conjg(u(2,c,l))*conjg(u(2,i,q))*J_in(b,d,m,k,l,q) *p                      &
  1261.  
  1262.                              -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  1263.  
  1264.                                 conjg(u(2,a,k))*u(2,c,l)*conjg(u(2,i,q))*K_in(b,d,m,k,l,q) *p                      &
  1265.  
  1266.                              -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*integral(k,q,l)* &
  1267.  
  1268.                                 u(2,a,k)*conjg(u(2,c,l))*u(2,i,q)*L_in(b,d,m,k,l,q)*conjg(p)  )
  1269.  
  1270.                  kn(k) = kn(k) + summ(k)
  1271.  
  1272.                   end if
  1273.  
  1274.                 end do
  1275.  
  1276.                end do
  1277.  
  1278.               end do
  1279.  
  1280.              end do
  1281.  
  1282.             end do
  1283.  
  1284.            end do
  1285.  
  1286.           end do
  1287.  
  1288.          end do
  1289.  
  1290.         end do
  1291.  
  1292.  
  1293.  
  1294.          Do k=1, estep
  1295.  
  1296.          kn(k) = kn(k) - n(k)/tau_x + pump_n(t,k)
  1297.  
  1298.          end do
  1299.  
  1300.  
  1301.  
  1302.      end subroutine n_deriv
  1303.  
  1304.  
  1305.  
  1306.  
  1307.  
  1308.  
  1309.  
  1310.  
  1311.  
  1312.      subroutine p_deriv(t, p, n, u, A_in, B_in, C_in, D_in, kp)
  1313.  
  1314.  
  1315.  
  1316.       implicit none
  1317.  
  1318.       complex(8), intent(in), dimension(2,2,2,estep,estep,estep) :: A_in, B_in, C_in, D_in
  1319.  
  1320.       real(8), intent(in) :: t
  1321.  
  1322.       complex(8), intent(in):: p
  1323.  
  1324.       complex(8), dimension(estep), intent(in) :: n
  1325.  
  1326.       complex(8), intent(in), dimension(2,2,estep) :: u
  1327.  
  1328.       complex(8), intent(out):: kp
  1329.  
  1330.       integer :: a,b,c,d,i,m,k,l,q
  1331.  
  1332.       complex(8):: sump
  1333.  
  1334.       real(8), external :: f0, integral , pump_p
  1335.  
  1336.  
  1337.  
  1338.  
  1339.  
  1340.       kp=(0.0,0.0)
  1341.  
  1342.       sump = (0.0,0.0)
  1343.  
  1344.  
  1345.  
  1346.         do k=1, estep
  1347.  
  1348.          do q=1, estep
  1349.  
  1350.           do l=1, estep
  1351.  
  1352.            do a=1,2
  1353.  
  1354.             do b=1,2
  1355.  
  1356.              do c=1,2
  1357.  
  1358.               do d=1,2
  1359.  
  1360.                do i=1,2
  1361.  
  1362.                 do m=1,2
  1363.  
  1364.                  
  1365.  
  1366.                  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
  1367.  
  1368.                   sump = (f0(k,q)*f0(k,q)*S**2/(2.0*pi**3))*                                                                &
  1369.  
  1370.                        ((sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(real(k)*pi/Lc)*integral(k,q,l)* &
  1371.  
  1372.                            u(1,a,k)*conjg(u(1,c,l))*conjg(u(1,i,q))*A_in(b,d,m,k,l,q)                        &
  1373.  
  1374.                        -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(real(k)*pi/Lc)*integral(k,q,l)* &
  1375.  
  1376.                            conjg(u(1,a,k))*u(1,c,l)*conjg(u(1,i,q))*B_in(b,d,m,k,l,q)                        &
  1377.  
  1378.                        +(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(real(k)*pi/Lc)*integral(k,q,l)* &
  1379.  
  1380.                            u(2,a,k)*conjg(u(2,c,l))*conjg(u(2,i,q))*C_in(b,d,m,k,l,q)                        &
  1381.  
  1382.                        -(sigma_z(a,b)*sigma_z(c,d)*sigma_z(i,m))*(real(k)*pi/Lc)*integral(k,q,l)* &
  1383.  
  1384.                            conjg(u(2,a,k))*u(2,c,l)*conjg(u(2,i,q))*D_in(b,d,m,k,l,q) )
  1385.  
  1386.                   kp = kp + sump
  1387.  
  1388.  
  1389.  
  1390.                   end if
  1391.  
  1392.                 end do
  1393.  
  1394.                end do
  1395.  
  1396.               end do
  1397.  
  1398.              end do
  1399.  
  1400.             end do
  1401.  
  1402.            end do
  1403.  
  1404.           end do
  1405.  
  1406.          end do
  1407.  
  1408.         end do
  1409.  
  1410.     kp= conjg(kp) - p/tau_c + pump_p(t)
  1411.  
  1412.    
  1413.  
  1414.     end subroutine p_deriv
  1415.  
  1416.   End module runge
  1417.  
  1418.  
  1419.  
  1420.  
  1421.  
  1422. !===================================================================================================
  1423.  
  1424.  
  1425.  
  1426. !Bat dau chuong trinh chinh.
  1427.  
  1428. !Chuong trinh nay se doc cac tham so dau vao: e(k) tu file. Doc cac gia tri t, h, f0 tu nguoi dung
  1429.  
  1430. ! Sau do se xuat ra file các gia tri cua p(t), n(k,t)
  1431.  
  1432.  
  1433.  
  1434.   Program process
  1435.  
  1436.     use global
  1437.  
  1438.     use runge
  1439.  
  1440.     implicit none
  1441.  
  1442.  
  1443.  
  1444.     complex(8) :: p                                
  1445.  
  1446.     complex(8):: p_out
  1447.  
  1448.  
  1449.  
  1450.     complex(8) , dimension(estep) :: n            
  1451.  
  1452.     complex(8), dimension(estep) :: n_out
  1453.  
  1454.  
  1455.  
  1456.     real(8) :: t  , t_max
  1457.  
  1458.     real(8) :: t_out                        
  1459.  
  1460.    
  1461.  
  1462.     complex(8),  dimension(2,2,estep):: u
  1463.  
  1464.     complex(8),  dimension(2,2,estep):: u_out  
  1465.  
  1466.  
  1467.  
  1468.     complex(8), dimension(2,2,2,estep,estep,estep) :: A_in, B_in, C_in, D_in, E_in, F_in,&
  1469.  
  1470.                                                                   G_in, H_in, I_in, J_in, K_in, L_in
  1471.  
  1472.     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,&
  1473.  
  1474.                                                  I_out, J_out, K_out, L_out
  1475.  
  1476.         real(8), external :: f0
  1477.  
  1478.     integer :: i
  1479.  
  1480.  
  1481.  
  1482.  sigma_z(1,1) = (1.0,0.0)
  1483.  
  1484.  sigma_z(1,2) = (0.0,0.0)
  1485.  
  1486.  sigma_z(2,1) = (0.0,0.0)
  1487.  
  1488.  sigma_z(2,2) =(-1.0,0.0)
  1489.  
  1490.  sigma_y(1,1) = (0.0,0.0)
  1491.  
  1492.  sigma_y(1,2) =-iu
  1493.  
  1494.  sigma_y(2,1) = iu
  1495.  
  1496.  sigma_y(2,2) = (0.0,0.0)
  1497.  
  1498.  
  1499.  
  1500.  
  1501.  
  1502. Print *,"Nhap t0"
  1503.  
  1504. read *, t
  1505.  
  1506.  
  1507.  
  1508. Print *, "Nhap h, voi h la buoc cua t"
  1509.  
  1510. read * , h
  1511.  
  1512.  
  1513.  
  1514. Print *, "Nhap t_max"
  1515.  
  1516. read *, t_max
  1517.  
  1518.  
  1519.  
  1520. Print *, "Nhap e(k)"
  1521.  
  1522.   open (unit =2, file="/home/nap/Desktop/list_E.txt")
  1523.  
  1524.   do i=1, estep
  1525.  
  1526.     read (2,*) e(i)
  1527.  
  1528.  
  1529.  
  1530.     print *, e(i)
  1531.  
  1532.   end do
  1533.  
  1534.  
  1535.  
  1536. p= (5.d-2,5.d-2)
  1537.  
  1538. n=(0.0,0.1)
  1539.  
  1540. u=(1.0,0.0)
  1541.  
  1542. A_in=(0.0,0.0)
  1543.  
  1544. B_in=(0.0,0.0)
  1545.  
  1546. C_in=(0.0,0.0)
  1547.  
  1548. D_in=(0.0,0.0)
  1549.  
  1550. E_in=(0.0,0.0)
  1551.  
  1552. F_in=(0.0,0.0)
  1553.  
  1554. G_in=(0.0,0.0)
  1555.  
  1556. H_in=(0.0,0.0)
  1557.  
  1558. I_in=(0.0,0.0)
  1559.  
  1560. J_in=(0.0,0.0)
  1561.  
  1562. K_in=(0.0,0.0)
  1563.  
  1564. L_in=(0.0,0.0)
  1565.  
  1566.  
  1567.  
  1568.  
  1569.  
  1570.  
  1571.  
  1572. open (unit=1, file="/home/nap/Desktop/result1.txt")
  1573. open (unit=3, file="result2.txt")
  1574. open (unit=4, file="result3.txt")
  1575.  
  1576.  
  1577.  
  1578. do while(t <= t_max)
  1579.  
  1580.  
  1581.  
  1582.  call rk4(t, h, e ,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,&
  1583.  
  1584.                               A_out, B_out, C_out, D_out, E_out, F_out,G_out, H_out, I_out, J_out, K_out, L_out, &
  1585.  
  1586.                                t_out, p_out, n_out, u_out)
  1587.  
  1588.   print *, t_out
  1589.  
  1590.  
  1591.  
  1592.   write (1,*)   n_out(2)
  1593.   write (3,*)   p_out
  1594.   write (4,*)   u(:,:,2)
  1595.  
  1596.   t = t_out
  1597.  
  1598.   p = p_out
  1599.  
  1600.   n = n_out
  1601.  
  1602.   u = u_out
  1603.  
  1604.   A_in = A_out
  1605.  
  1606.   B_in = B_out
  1607.  
  1608.   C_in = C_out
  1609.  
  1610.   D_in = D_out
  1611.  
  1612.   E_in = E_out
  1613.  
  1614.   F_in = F_out
  1615.  
  1616.   G_in = G_out
  1617.  
  1618.   H_in = H_out
  1619.  
  1620.   I_in = I_out
  1621.  
  1622.   J_in = J_out
  1623.  
  1624.   K_in = K_out
  1625.  
  1626.   L_in = L_out
  1627.  
  1628.  
  1629.  
  1630. end do
  1631.  
  1632.  
  1633.  
  1634.  
  1635.  
  1636. end program process
  1637.  
  1638.  
  1639.  
  1640.  
  1641.  
  1642.  
  1643.  
  1644. !=================================================================================  
  1645.  
  1646.  function f0(a,b)
  1647.  
  1648.   use global
  1649.  
  1650.   implicit none
  1651.  
  1652.    real(8) :: f0, k, q, l, Ec_k, Ex_k, Ec_q, Ex_q, Ec_l, Ex_l, Ec_0, Ex_0
  1653.  
  1654.    integer :: a,b
  1655.  
  1656.  
  1657.  
  1658.       k = real(a)*pi/Lc
  1659.  
  1660.       q = real(b)*pi/Lc
  1661.  
  1662.       l=k-q
  1663.  
  1664.                                                          
  1665.  
  1666.  
  1667.  
  1668.    Ec_k = sqrt(Ec0**2+ (hbar*k*k/sqrt(13.5))**2)
  1669.  
  1670.  
  1671.  
  1672.    Ex_k = E0 + hbar*(k**2)/(2.0*(m_e + m_h))
  1673.  
  1674.  
  1675.  
  1676.    Ec_q = sqrt(Ec0**2+ (hbar*q*q/sqrt(13.5))**2)
  1677.  
  1678.  
  1679.  
  1680.    Ex_q = E0 + hbar*(q**2)/(2.0*(m_e + m_h))
  1681.  
  1682.  
  1683.  
  1684.    Ec_l = sqrt(Ec0**2+ (hbar*l*l/sqrt(13.5))**2)
  1685.  
  1686.  
  1687.  
  1688.    Ex_l = E0 + hbar*(l**2)/(2.0*(m_e + m_h))
  1689.  
  1690.  
  1691.  
  1692.    Ec_0 = sqrt(Ec0**2)
  1693.  
  1694.  
  1695.  
  1696.    Ex_0 = E0
  1697.  
  1698.  
  1699.  
  1700.    f0 = 6.0*10.0*100.0 * 1.0/ ( sqrt(1.0+(hbar_omega/(2.0*hbar*(Ec_k-Ex_k)))**2)*  &
  1701.  
  1702.                              sqrt(1.0+(hbar_omega/(2.0*hbar*(Ec_q-Ex_q)))**2)* &
  1703.  
  1704.                              sqrt(1.0+(hbar_omega/(2.0*hbar*(Ec_l-Ex_l)))**2)* &
  1705.  
  1706.                              sqrt(1.0+(hbar_omega/(2.0*hbar*(Ec_0-Ex_0)))**2)*S )
  1707.  
  1708.  
  1709.  
  1710. End function f0
  1711.  
  1712.  
  1713.  
  1714. !=======================================================================
  1715.  
  1716. ! ham nay rut gon bieu thuc tinh k, p, l
  1717.  
  1718.  
  1719.  
  1720. Function integral(a,b,d)  
  1721.  
  1722.  use global
  1723.  
  1724.  implicit none
  1725.  
  1726.  real(8) :: integral, k, q, l
  1727.  
  1728.  integer :: a,b,d
  1729.  
  1730.  
  1731.  
  1732.       k = real(a)*pi/Lc
  1733.  
  1734.       q = real(b)*pi/Lc
  1735.  
  1736.       l = real(d)*pi/Lc
  1737.  
  1738.  
  1739.  
  1740.     integral = q*l/(sqrt(4.0*k*k*q*q-(k*k+q*q-l*l)**2))
  1741.  
  1742.  
  1743.  
  1744. end function integral
  1745.  
  1746.  
  1747.  
  1748.  
  1749.  
  1750.  
  1751.  
  1752.  
  1753.  
  1754. ! ======================================================================
  1755.  
  1756. ! Ham bom cua p(t)
  1757.  
  1758.  
  1759.  
  1760. Function pump_p(t)
  1761.  
  1762.  use global
  1763.  
  1764.  implicit none
  1765.  
  1766.  
  1767.  
  1768.  real, parameter :: delta_P = 5.d-2  , t_p= 3.0    ! gia tri tam thoi de test
  1769.  
  1770.  real(8)::     pump_p, t
  1771.  
  1772.  
  1773.  
  1774.  pump_p = delta_P*exp(-2.0*log((2.0*t**2+0.01)/t_p**2))
  1775.  
  1776.  
  1777.  
  1778. end function pump_p
  1779.  
  1780.  
  1781.  
  1782. !======================================================================
  1783.  
  1784. ! Ham bom cua n(t)
  1785.  
  1786.  
  1787.  
  1788. function pump_n(t, a)
  1789.  
  1790.  
  1791.  
  1792.  use global
  1793.  
  1794.  implicit none
  1795.  
  1796.  
  1797.  
  1798.  real, parameter:: P_0 =5.d-2 , t_p =3.0 , delta_E =0.25    ! gia tri tam thoi
  1799.  
  1800.  real(8) :: pump_n , t ,k
  1801.  
  1802.  integer :: a
  1803.  
  1804.  
  1805.  
  1806.  k = real(a)*pi/Lc
  1807.  
  1808.  pump_n = P_0 *((E(a+1)-1.61579d3)/delta_E)*    exp(-2.0*log((2.0*t**2+0.01)/t_p**2))
  1809.  
  1810.  
  1811.  
  1812.  end function pump_n
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!