function[E_value, JFN_value exit_flag] = Efield2(guess,psi,beta,d,p,T,gamma_SE,gas_type)
     


%A_FN = 6.2e-6;% (L)
%B_FN  = 6.85e9;% (L)
A_FN = 1.5414*10^-6; %(V.A)
B_FN = 6.8308*10^9; %(V.A)
t_2_y = 1.1;
k= 1.3806488*10^-23;
e =1.6021765*10^-19;
epsilon = 8.8541878 *10^-12;
AMU = 1.6605e-27;
sigma_CE =9.67e-19;
C_FN = A_FN*beta^2/(psi*t_2_y)*exp((3.79*10^-5)^2*B_FN/sqrt(psi));
D_FN = 0.95*B_FN*psi^(3/2)/beta;
switch gas_type
    
    case 1 
        return
    case 2 % Air
        A_p = 15; B_p = 365;
        m = AMU*28.96;
        IP = 14.9;
        p_star = D_FN*B_p^(-2);
        L = 1/(p_star*A_p);
    case 3 % Argon
        C_p = 29.2*(10^2/133.32); D_p = 26.6*(10^2/133.32)^(1/2);  
        ap_macro  = @(E) p*C_p*exp(-D_p*sqrt(p/E));
        m=39.978*AMU;
        IP = 15.76;
        p_star = D_FN*D_p^(-2);
        L = 1/(p_star*C_p);
    case 4 % Helium  % problem at 535 start_dis = 7
        C_p = 4.4*(10^2/133.32); D_p = 14 * (10^2/133.32)^(1/2);
        ap_macro  = @(E)p* C_p*exp(-D_p*sqrt(p/E));
        m= 4.002*AMU;
        IP = 24.587;
        p_star = D_FN*D_p^(-2);
        L = 1/(p_star*C_p);
    case 5 % Nitrogen
        A_p = 12 ; B_p = 342;
        m = 28.01*AMU;
        IP = 14.53;
        p_star = D_FN*B_p^(-2);
        L = 1/(p_star*A_p);
    case 6 % Oxygen
        A_p = 14; B_p = 341;
        m =31.99*AMU;
        IP =13.61;
        p_star = D_FN*B_p^(-2);
        L = 1/(p_star*A_p);
    case 7 % Xenon
        C_p = 65.3*(10^2/133.32); D_p = 36.1*(10^2/133.32)^(1/2);   
        ap_macro  = @(E) p*C_p*exp(-D_p*sqrt(p/E));
        m =131.29*AMU;
        IP = 12.13;
        p_star = D_FN*D_p^(-2);
        L = 1/(p_star*C_p);
end

%Air Correction
if gas_type ~= 3 && gas_type ~= 4 && gas_type ~= 7
ap_macro = @(E) A_p*p*10^2/133.33*exp(-B_p*p*10^2/(133.3*E));
end

alpha  = @(E) ap_macro(E)*(1 - exp(-((E*d/IP -1)/3.1)^0.8));

JFN =@(E) C_FN*E^2*exp(-D_FN/E);


T_star = ( (pi*p_star*m*sigma_CE)/( D_FN*2*e*k ) )*...
           ( ((D_FN*C_FN*L)/(2*epsilon))^2 ) ;
J_0 = C_FN*D_FN^2;
   
E_n = @(E) E/D_FN;
JFN_n = @(E) JFN(E)/J_0;

T_n = T/T_star;
alpha_n = @(E) alpha(E)*L;
p_n = p/p_star;

d_n = d/L;    
pd = p_n*d_n;
x_0 = @(E)( sqrt(1+8*E_n(E))-1 )/( 4*E_n(E) );


F_br =@(E) sqrt( (T_n*E_n(E))/(p_n*(d_n^2)) ) *...
    ( exp(1/E_n(E))*(1-gamma_SE*(exp(alpha_n(E)*d_n)-1)) / ...
    (exp(alpha_n(E)*d_n) - 1) );

G =@(E) F_br(E) - exp(x_0(E))*(1+2*E_n(E)*x_0(E))/x_0(E);


[E_value, ~, exit_flag] = fzero(G,guess);
JFN_value = JFN(E_value);