rho1=1;
    h=1;
    sigma1=sqrt(10^-3);
    alpha=2;
    beta=(2500);
    d=40; % assumption
    r=[50 ,55,60,65,70];
    n=100; % number of  data points 
    x_axis=linspace(0,35,8)
    Pt_values=(10^-3)*10.^(x_axis/10);
    %n - no of values for z 
    z0=1; % cdf meh z = 1 daal rahe hai
    startpow=10;
    endpow=1000;
    noofpow=30;
    out1=zeros(1,8);
    k=10^3.7;
    cindex = 1;
    for index =1:8
        Pt=Pt_values(index);
        c= (rho1 * Pt)/sigma1^2;  
        sum=0;
        for i=1:5
            r0=sqrt(r(i)^2+d^2);            
            z=z0;
            %const j
     
            %(exp(-((r0^alpha)*z0)/(c*beta))/(r0^alpha))
            %
            j=(1/(r0^alpha));

            sum=sum +((exp(-1*(((r0^alpha)*z)/(c *beta))))/(r0^alpha) - j) ;
        end
        index;
        out1(index)=sum*(-0.2)*k ;   
    end
    out1
     
    
     
    plot(x_axis,out1);
    semilogy(x_axis,out1);
    ylabel("Outage Probability")
    xlabel("c in decibel, where c is gamma1 dash")
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    r=[50 ,55 ,60 ,65 ,70];
    alpha1=3;
    alpha2=4;
    lambda1=1.12201
    lambda2=1.12201
    omega1=1
    omega2=1
    q=1;
    % Calculations for SNR2 |
    %Expected value of higher order moments of h.g where h and g are rician
    %random variables 
    %qth order moment
     
    alpha11=2.33
    alpha12=1.4
     
    M=mu_k(1,lambda1,lambda2,omega1,omega2,alpha1,alpha2);
     
    %defining phis
    no_of_irs_ele=4
    % lets call 2 paramters , alpha 11 , alpha 12 
     
    phi2=mu_k(2,lambda1,lambda2,omega1,omega2,alpha11,alpha12) /mu_k(1,lambda1,lambda2,omega1,omega2,alpha11,alpha12)
    phi3=mu_k(3,lambda1,lambda2,omega1,omega2,alpha11,alpha12) /mu_k(2,lambda1,lambda2,omega1,omega2,alpha11,alpha12)
    phi4=mu_k(4,lambda1,lambda2,omega1,omega2,alpha11,alpha12) /mu_k(3,lambda1,lambda2,omega1,omega2,alpha11,alpha12)
    mu_1=no_of_irs_ele*mu_k(1,lambda1,lambda2,omega1,omega2,alpha11,alpha12)
    %defining 'a ' constants 
    a3 = ( 4*phi4 - 9*phi3 +6*phi2 - mu_1) /(-phi4 + 3*phi3 - 3*phi2 + mu_1 )
    a2 = ((a3/2)  * ( phi4 - 2*phi3 + phi2) ) + 2*phi4 - 3*phi3 + phi2 
    %a6 = ((a3*(phi2 - mu_1) + 2*phi2 - mu_1 )/a2 ) -3
    a6 = (1/a2)*(a3*(phi2 - mu_1) + 2*phi2 - mu_1) -3
    a7 = sqrt((a6+2)^2 - ((4*mu_1*(a3+1))/a2))
    a5= (a6-a7)/2
    a4 = ( a6+ a7 )/2
    a1 = gamma(abs(a3+1)) / ( a2 * gamma(abs(a4+1))*gamma(abs(a5+1)))
     
     
    %defining the CDF function of meijer G ; consider it to be 5 different
    %independent distributions and act accordingly 
    % make it into a summation format 
     
    a_vec=[1];
    b_vec=[a3+1];
    c_vec=[a4+1,a5+1]
    d_vec=[0];
    dist2=(500)^alpha2;
     
    sigma2=((10)^(-7.5));
    no_of_points=8;
    maxpow=35;
    minpow=0;
    out2=zeros(1,no_of_points);
    %x_axis=linspace(minpow,maxpow,no_of_points);
    Pt2_values=10.^(x_axis/10)*(10^-3);
    rho2=1;
     
     
     
     
    for index=1:no_of_points
        curr_pow=Pt2_values(index);
        gamma2=curr_pow/(sigma2^2);
        sum=0;
        for i=1:5
            dist1=r(i)^alpha1
            dist2
            rho2*gamma2
            z=(1/a2)*sqrt((dist1*dist2)/(rho2*gamma2))
            % ai ,a2 * value 
            value=a1*a2*meijerG(a_vec,b_vec,c_vec,d_vec,z);
            sum=sum + real(value);
        end
        %(10^(-(index-1)/3))*
        out2(index)=0.2*sum;
    end
     
    out2
    plot(x_axis,out2)
    title('outage prob 2 ')
    semilogy(x_axis,out2)
    title('outage prob 2 semilog')
     
    out4=1.-out1;
    out1
    out2
    out10=out1+out2
    out11=out1.*out2
    out3=out1+out2-(out1.*out2)
    plot(x_axis,out3)
    title('final outage  ')
     
    semilogy(x_axis,out3)
    title('final outage semilog ')
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     
     

    function expected_moment = mu_k(q,lambda1,lambda2,omega1,omega2,alpha11,alpha12)
        e1=laguerreL(q,0,- (((lambda1)^2)/omega1));
        e2=laguerreL(q,0,-(((lambda2)^2)/omega2));
        Exp =((factorial(q))^2 /  (alpha11*alpha12)^q)*e1*e2;
        expected_moment = Exp; % Mean value of the elements in q
    end