Indexing inside double for loops

guilherme stahlberg

2019 8 5

0 个回答

9 次查看（30 天）

0 个投票

Hey everyone,

I'm trying to index a vector inside a double for loop, but can't figure out what's going wrong. As one can see from the code below, I'm using the En function that depends on ft and s to calculate U and a0 which are used to calculate En. At the next step, En is changing so U and a0 will also change, but I can't figure how to index them correcty inside their loops. Can anyone help me sort this out? Thanks in advance.

Kind regards

format long
clc
%% Variables
    fs=1/10; 
    s=0:fs:1;
    ft = -2980:fs:2980;
    alpha = 0.01;
    ai = 0.4;
    ai2 = 0.6;
    af = 1;
    w0 = 0.02;
    mi = 6;
    mmi=[0 mi;mi 0];
    E1 = 0.03;
    E2 = 0.01;
    psii=[sqrt(ai);sqrt(ai2)];
        
%% Functions
    g = 1./(1+exp(-alpha.*ft));
        
    f = ai*(1-g)+af*g;
        
%% Electric Field
    E0 = alpha.*(af-ai).*exp(alpha.*ft).*sin(w0.*ft)./(mi.*(1+exp(alpha.*ft)).*sqrt((1-ai+(1-af).*exp(alpha.*ft)).*(ai+af.*exp(alpha.*ft))));
                
%% New Electric Field
for k=1:length(s)
        
        if k == 1
            
            En(k)=E0(k);
            
        end
        
    % Time Operador
        
    H0=[-1i*E1*fs/2 0;0 -1i*E2*fs/2];
    H01=expm(H0);
    Z=[0 1;1 0];
    [V,D]=eig(Z);
        
    for i = 1:length(ft)
    
        H(:,:,i)=1i*mi*En(i)*fs*D;
        H1(:,:,i)=expm(H(:,:,i));
        Hint(:,:,i)=V*H1(:,:,i)*transpose(V); 
        U0(:,:,i)=H01*Hint(:,:,i)*H01;
        
        if i==1
   
            U(:,:,i)=U0(:,:,i);
        
            else
            
            U(:,:,i)=U0(:,:,i)*U(:,:,i-1);
            
        end
            
    psi(:,i)=U(:,:,i)*psii;
    psi1(i)=psi(1,i);
    psi2(i)=psi(2,i);
            
    end
    psif=psi(:,59601);
    Uf1=U(:,:,59601);
  
    % a0
    for j = 1:length(ft)
   
        Ui(:,:,j)=U(:,:,j);
        Uf(:,:,j)=Uf1*ctranspose(Ui(:,:,j));    
        brapsii=ctranspose(psii);
        brapsif=ctranspose(psif);
        Udaga=ctranspose(Uf1);
        a0(j)=-2*imag(brapsii*Udaga*psif*brapsif*Uf(:,:,j)*mmi*Ui(:,:,j)*psii);
    end
    % New Electric Field
    
    S=1;
    a0int=S.*a0.^2;
    gamao(k)=trapz(ft,a0int);
    fe(k)=En(k)/(600*fs);
    gint=S.*a0.*fe;
    gaminha(k)=trapz(ft,gint);
        
    En(k+1)=En(k)+(fe(k)-gaminha(k)*a0(k)/gamao(k))*fs;
        
    figure
    plot(ft,En)
    grid on 
    
end