0
$\begingroup$

I want to plot sigma according to lambda.Sigma has to generate from loop code. but I can't extract these data. please help me with this problem. I need to calculate the sigma value for a given range of lambda, for which I have to define the masculine epsilon and alpha and W matrices. Finally, I can use this information to get sigma. Then plot a sigma for all the loops in a curve. At first, I was import data from the file and introduced some permanent:

Clear["Global`*"];
        SetDirectory[NotebookDirectory[]];
        NN = Import["n-Au.txt", "Table"];
        KK = Import["k-Au.txt", "Table"];
        NN11 = NN[[All, 1]]*(10^(-4)); NN22 = NN[[All, 2]];    (*n*)
        KK22 = KK[[All, 2]];
        N1 = Table[{NN11[[i]], NN22[[i]]}, {i, 1, Length[NN11]}];
        K1 = Table[{NN11[[i]], KK22[[i]]}, {i, 1, Length[NN11]}];
        N2 = Interpolation[N1];
        K2 = Interpolation[K1];
        c = 2.99*10^10;
        a = 10.0*10^(-7);
        bx = 0.0;
        by = 0.0;
        bz = 8000;
        Np = 2;
        A = 3/4.0;
        omegap = 1.386*10^16;
        gammabulk = 4.584*10^13;
        vf = 1.39*10^8;
        gamma = gammabulk + (A*vf)/a;
        epsilonm = 1.7;
        e = 4.8032*10^(-10);
        me = 9.1*10^(-28);
        B = bz;
        omegac = (e*B)/(me*c);
        phi = 0.0;
        theta = Pi/4.0;
        phi0 = 0.0;
        theta0 = 0.0;
        deltalambda = 0.01*10^(-4);
        ii = 1;
        r[i_, j_] := r[i] - r[j];
        r[1] = {0.0, 0.0, 0.0};
        r[2] = {3.0*a, 0.0, 0.0};

here I have introduced some matrix that I need in my calculation

        For[lambda = 0.2*10^(-4), lambda <= 0.8*10^(-4), 
         lambda = lambda + deltalambda,
         omega = 2.0*Pi*c/lambda;
         k = omega*Sqrt[epsilonm]/c;
         kvector = k {Sin[theta]*Cos[phi], Sin[theta]*Sin[phi], Cos[theta]};
         realepsilonbulk = N2[lambda]^2 - K2[lambda]^2;
         imaginaryepsilonbulk = 2*N2[lambda]*K2[lambda];
         G = omegap^2/(omega^2 + I*gamma*omega)^2;
         epsilonxy = I*G*omegac*bz;
         epsilonxz = -I*G*omegac*by;
         epsilonyx = -I*G*omegac*bz;
         epsilonyz = I*G*omegac*bx;
         epsilonzx = I*G*omegac*by;
         epsilonzy = -I*G*omegac*bx;
         reepsilonxx = 
          realepsilonbulk + 
           Re[omegap^2/(omega^2 + I*omega*gammabulk) - omegap^2/(
             omega^2 + I*omega*gamma)];
         imepsilonxx = 
          imaginaryepsilonbulk + 
           Im[omegap^2/(omega^2 + I*omega*gammabulk) - omegap^2/(
             omega^2 + I*omega*gamma)];
         epsilonxx = reepsilonxx + I*imepsilonxx;
         epsilon = {{epsilonxx, epsilonxy, epsilonxz}, {epsilonyx, epsilonxx, 
            epsilonyz}, {epsilonzx, epsilonzy, epsilonxx}};
         F = -3*I*a^3*omegap^2/(omega*(omega + I gamma)^2)*
           epsilonm/(epsilonxx + 2*epsilonm)^2;
         alphaxx = a^3*(epsilonxx - epsilonm)/(epsilonxx + 2*epsilonm);
         alpha = {{alphaxx, -F*omegac*bz, F*omegac*by}, {F*omegac*bz, 
            alphaxx, -F*omegac*bx}, {-F*omegac*by, F*omegac*bx, alphaxx}};
         thetahatr = {{(1/Sqrt[2])*(Cos[theta]*Cos[phi] - I*Sin[phi])}, {(1/
               Sqrt[2])*(Cos[theta]*Sin[phi] + I*Cos[phi])}, {(1/
               Sqrt[2])*(-Sin[theta])}};
         thetahatl = {{(1/Sqrt[2])*(Cos[theta]*Cos[phi] + I*Sin[phi])}, {(1/
               Sqrt[2])*(Cos[theta]*Sin[phi] - I*Cos[phi])}, {(1/
               Sqrt[2])*(-Sin[theta])}};
         E0 = 1.0;
         El1r = (E0*(Exp[I*(kvector.r[1])])) thetahatr;
         El2r = (E0*(Exp[I*(kvector.r[2])])) thetahatr;
         El1l = (E0*(Exp[I*(kvector.r[1])])) thetahatl;
         El2l = (E0*(Exp[I*(kvector.r[2])])) thetahatl;
         For[i = 1, i <= Np, i++,
          For[j = 1, j <= Np, j++,
           For[be = 1, be <= 3, be++,
            For[al = 1, al <= 3, al++,
             If[i == j, wdiagonal = Table[0, {al, 1, 3}, {be, 1, 3}],
              If[i == 1,
        w12 = Table[
          E^(I*k*Norm[
               r[i] - r[
                 j]])*((k^2*(KroneckerDelta[al, be]*Norm[r[i] - r[j]]^2 - 
                  r[i, j][[al]]*r[i, j][[be]]))/(epsilonm*
                Norm[r[i] - r[j]]^3) + 
                  ((1 - 
                  I*k*Norm[r[i] - r[j]])*(3*(r[i, j][[al]]*r[i, j][[be]]) - 
                  KroneckerDelta[al, be]*Norm[r[i] - r[j]]^2))/(epsilonm*
                Norm[r[i] - r[j]]^5)), {al, 1, 3}, {be, 1, 3}]
                   , w21 = Table[((k^2*(Norm[r[i] - r[j]]^2*KroneckerDelta[al, be] - 
                  r[i, j][[al]]*r[i, j][[be]]))/(epsilonm*
                Norm[r[i] - r[j]]^3) + 
                  ((1 - 
                  I*k*Norm[r[i] - r[j]])*(3*(r[i, j][[al]]*r[i, j][[be]]) - 
                  Norm[r[i] - r[j]]^2*KroneckerDelta[al, be]))/(epsilonm*
                Norm[r[i] - r[j]]^5))*Exp[I*k*Norm[r[i] - r[j]]], 
             {al, 1, 3}, {be, 1, 3}] 
                   ]]]]]];

then calculated some value that I need:

         wmatrixblock = ArrayFlatten[{{wdiagonal, w12}, {w21, wdiagonal}}];
         wmatrixblockdotalpha = 
          ArrayFlatten[{{wdiagonal, alpha.w12}, {alpha.w21, wdiagonal}}];
         Amatrix = IdentityMatrix[6] - wmatrixblockdotalpha;
         alphae1r = alpha.El1r;
         alphae2r = alpha.El2r;
         alphaer = ArrayFlatten[{{alphae1r}, {alphae2r}}];
         alphae1l = alpha.El1l;
         alphae2l = alpha.El2l;
         alphael = ArrayFlatten[{{alphae1l}, {alphae2l}}];
         dblockr = {{d11r}, {d12r}, {d13r}, {d21r}, {d22r}, {d23r}};
         Amatrix.dblockr == Transpose[alphaer];
         solvr = LinearSolve[Amatrix, alphaer];
         dblockl = {{d11l}, {d12l}, {d13l}, {d21l}, {d22l}, {d23l}};
         Amatrix.dblockl == Transpose[alphael];
         solvl = LinearSolve[Amatrix, alphael];
         d1r = {solvr[[1]], solvr[[2]], solvr[[3]]};
         d2r = {solvr[[4]], solvr[[5]], solvr[[6]]};
         d1l = {solvl[[1]], solvl[[2]], solvl[[3]]};
         d2l = {solvl[[4]], solvl[[5]], solvl[[6]]};
         sigmar = ((4*Pi*omega)/((E0^2)*Sqrt[epsilonm]*
               c))*(Im[(Transpose[(Inverse[Conjugate[alpha]].Conjugate[
                    d1r])].d1r + 
               Transpose[(Inverse[Conjugate[alpha]].Conjugate[d2r])].d2r)]);
         sigmal = ((4*Pi*omega)/((E0^2)*Sqrt[epsilonm]*
               c))*(Im[(Transpose[(Inverse[Conjugate[alpha]].Conjugate[
                    d1l])].d1l + 
               Transpose[(Inverse[Conjugate[alpha]].Conjugate[d2l])].d2l)]);
         sigma = (sigmal - sigmar)/(Pi*(a^2));
         ]
$\endgroup$
  • 1
    $\begingroup$ Don't use For. Better use Do. But for generating a table, use, well, Table. $\endgroup$ – Henrik Schumacher Dec 1 '18 at 12:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.