API tools faq
paste
Advertisement
Guest User

Untitled

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