% hybrid 3D (B in XY plane)

% +---------+
% FORESHOCK.M
% +---------+
%
% Original source code written 16/11/1999 by Stuart Knock in Fortran90s
%
% Calculates relaxed velocity distribution of a given location in
% the foreshock region of a planetary/interplanetary shock.
% For details of equations used in calculating cut off distribution
% see I.Cairns(1987) "The Electron Distribution Function Upstream
% From the Earth's Bow Shock"
%
% Modification 22/11/99 (working)
% 23/11/99 added calculate_emissivities subroutine (working)
% 25/11/99 modified to calculate over array of positions(working)
%
% 22/1/00 Incorporated shock drift acceleration as the mechanism of
% electron acceleration and therefore its associated distribution 
% function produced at the shock. 
%
% 11/2/00 Changed number density for beam to include only those 
% particles  actually in the beam  ie. without the background 
% particles incorporated. 
%
% +-+-+ Modifications from this point made by Zdenka Kuncic +-+-+
%
% July 2000: optimization begun; several modifications to code to
% improve run-time performance, followed by compilation options.
%
% August 2000: modifications made to apply code to model Earth's
% foreshock (program renamed accordingly as foreshock.f90);
% analytical model for curvature constant implemented.
%
% Sept. 2000: self-consistency improved by calculating MHD shock
% jump conditions (Iver's routine cgroots used) and by adding a
% cross-shock potential model derived analytically using the
% electron momentum equation and the shock jump conditions.
%
% March 2001: replaced Maxwellian background (and incorrect
% reduced loss-cone distribution) with a more realistic Gaussian
% core + kappa halo distribution.
%
% August 2001: extended calculations over entire shockfront, ie
% include flipside of foreshock where R < 0 and thetabn > 90.
%
% December 2001: modifications to improve code performance in
% producing output data taking advantage of f90 array procedures.
%
% Sept/Oct 2002: refinement of advective free energy density flux
% calculation used in emissivities routine to include factor 
% deltavb/vb and contribution from vd advection in x-direction.
% Addition of Langmuir wave electric field strength calculation.
%
% May 2003: modified Langmuir wave damping rate to include ES decay.
%
% October 2003: minor modifications to input data to generalize to
% planetary foreshocks. 
%
% +-+-+- Modifications from this point made by Patrick Oppel +-+-+-
%
% Jan/Feb 2020: translating into Matlab code
% February 2020: added linear mode conversion efficiency and resulting
% emissivity, added total emitted power and observed power levels
%
% Mar/Apr 2020: adoption to bolide foreshock

% clearvars -except Rnxn f_hybrid a_hybrid b_hybrid n_hybrid f_z_og f_z_og_flat f_z_h f_z_h_flat vel_h vc_h vc_z vel
% clearvars -except RnxnF Rnxn s
% clearvars -except RpnF RpF s
% clearvars -except RxRangeF RxRange sRx
% clearvars -except NeF Nen sNe
clear all
%% Global variables
init

% Ne = Nen;
% Ne = 1e11;
% Rrange = RxRange(1)*Rp;
% xrange = RxRange(2)*Rp;

% Rp = RpF;
% Rrange = 200*Rp;   % spatial extent of source region
% xrange = 10*Rp;

% RR = Rnxn(1);
% xR = Rnxn(2);

% Physical constants
c = 2.9979e8;           % speed of light [m/s]
mi = 1.672614e-27;      % ion mass (assuming protons) [kg]
me = 9.109558e-31;      % electron mass [kg]
kbolt = 1.380622e-23;   % Boltzman constant [m^2 kg s^-2 K^-1]
e = 1.602e-19;          % electron charge [C]
mu0 = 1.25664e-6;       % permeability of free space [H m^-1]

if Matot < 1
    mi = 31*mi;
end

% dependent physical variables
ve = sqrt(kbolt*Te/me);     % electron thermal velocity
vec3 = (ve/c)^3;
vd = -vsw*sind(theta);     % ExB drift velocity
% vd = -vsw*sind(theta);     % ExB drift velocity
ctheta = cosd(theta);
stheta = sind(theta);
e21 = ecc*ecc-1;    
vmax = c;              % max velocity to which distributions are calculated
vbmin = 5*ve;          % minimum beam velocity      
% vbmin = ve;          % minimum beam velocity      

%% Program

Uc = 2;
cs = sqrt(5/3*kbolt*(Te+Ti)/mi);    % sound speed
pram = Ne*mi*vsw^2;                 % solar wind ram pressure
Ma = Matot;

Ms = vsw*abs(stheta)/cs;
% fcnX = cgroots2(cs,theta,Matot);
gam = 5/3;

if Ma > 1 
    na = -((gam-1)/(gam+1) + (gam+(gam+2))./((gam+1)*Ma.^2) + 2./((gam+1)*Ms.^2));
    nb = (gam*(1+1) - 2 + ((gam+1)./Ma.^2 + 4./Ms.^2))./((gam+1)*Ma.^2);
    nc = -1/((gam+1)*Ma.^4).*(gam - 1 + 2/Ms.^2);

    np = (3*nb - na.^2)/3;
    nq = (2*na.^3)/27 - na.*nb/3 + nc;
    nD = (np/3).^3 + (nq/2).^2;                                        

    nu = (-nq/2 + sqrt(nD)).^(1/3);
    nv = (-nq/2 - sqrt(nD)).^(1/3);                            
    
    % density jump across shock
    fcnX = nu + nv - na/3;
    as = 10.3*Rp*(3.4*fcnX+0.4)*(pram/1.8e-9)^(-1/6);
%     fcnX = cgroots(cs,theta,Matot);
else
    % normal gas dynamic shock
%     fcnX = (gam-1)/(gam+1)+2/(gam+1)/Ms^2;
%     fcnX = 1./((gam-1)/(gam+1)+2/(gam+1)/Ms^2);
%     fcnX_b = 1/(2*gam/(gam+1)*Ms^2);
%     pressure jump
%     fcnX = 1/(2*gam/(gam+1)*Ms^2);

%     density jump
%     fcnX = ((gam-1)*Ms^2+2)/((gam+1)*Ms^2);
    fcnX = (1+5*Ms^(-2))/6;
    as = 10.3*Rp*(3.4*fcnX+0.4)*(pram/1.8e-9)^(-1/6)+1;
end

% meteor model (Anderson, 2006)
% as = (fcnX/(1+sqrt(2*fcnX))+1)*Rp;
% check curvature constant without Matot
% curvature constant of bow shock model
% as = 10.3*Rp*(3.4*fcnX+0.4)*(pram/1.8e-9)^(-1/6);
% as_b = 10.3*Rp*(3.4*fcnX_b+0.4)*(pram/1.8e-9)^(-1/6);
% as = 10.3*Rp*(3.4*fcnX_b+0.4)*(pram/1.8e-9)^(-1/6);
% b = (as/(24.6+76.5/Matot^2)^2)*(pram/1.8e-9)^(1/3)/Rp^2;
b = as/(24.6)^2*(pram/1.8e-9)^(1/3)/Rp^2;
% b_delta = as_delta/(24.6)^2*(pram/1.8e-9)^(1/3)/Rp^2;
% b = as_b/(24.6)^2*(pram/1.8e-9)^(1/3)/Rp^2;
% b_deltaP = as_delta/100*(pram/1.8e-9)^(1/3)/Rp^2;

% tangent point of shock and B-field
if ecc == 1
    Yt = ctheta/2/b/stheta;
    Xt = as - b*Yt^2;
else
    cot2theta = (ctheta/stheta)^2;
    dum = e21 + 1/(as*b);
    aa = 1 - cot2theta*e21;
    bb = as*(cot2theta - 1/e21)*dum;
    nc = as/b/e21*(1-as*b/4*cot2theta*dum^2);
    Xt = (-bb-sqrt(bb^2-4*aa*nc))/(2*aa);
    Yt = ctheta/stheta*(Xt*e21-as/2*dum);
end

%% shockfront
xres = 8e3;

dR = 2*Rrange/(RR-1);
dx = xrange/xR;

% parameters for efficiency calculation
coneF = 0.25*pi; %solid angle of emission
coneH = 2*pi;
gamma_t = 1+3*Ti/Te;
vcrit = 1.5*sqrt(mi/(gamma_t*me))*ve;
temp4gammaLS = 8/3*vcrit*ve;
temp4zetaF = 0.2*(ve/vcrit)^2;
sq3 = sqrt(3);
erftemp1 = (ve*sq3/c)+(ve/vcrit);
erftemp2 = (ve*sq3/c)-(ve/vcrit);
temp4phiF = 72*exp(-Uc^2)/(Uc*sqrt(pi))*sq3*vec3;
temp4phiH = 12/5*sq3/gamma_t*vcrit/ve*vec3;


%% get_results

%% background distribution
vcore = sqrt(2)*ve;
vhalo = sqrt(kappa)*vcore;
ncore = Ne/(1+nhnc);
nhalo = nhnc*ncore;

%%
% vcRx = zeros(RR,xR);
vbRx = zeros(RR,xR);
% vbRx2 = zeros(RR,xR);
NbRx = zeros(RR,xR);
dvbRx = zeros(RR,xR);
% dvbRx2 = zeros(RR,xR);
Wfr = zeros(RR,xR);
l = zeros(RR,xR);
% thetabn = zeros(RR,xR,11,1001);
% 250
% B2B1 =  3.512327529379490*ones(1,1001);
% 249
% B2B1 =  3.512357386071995*ones(1,1001);
B2B1step = false;
        
tic
dZ = 1*Rp;
FF_nlt = 0;

loop = true;
Zn = 0;

M_f_constraint = 0;

iter0 = 10^(floor(log10(abs(Rp))));
itermin = iter0/1e5;

while loop
    Zn = Zn+1;
    if Zn == 1
        Z = 0;
    else
        Z = dZ*(Zn-1);
    end
    % spatial coordinate parallel to B
    for Rn = 1:RR
%     for Rn = 60
%     for Rn = Rnxn(1)
%         fprintf('Rn: %d of %d \n', Rn, RR)
        R = -Rrange + (Rn-1)*dR;                 

    % spatial coordinate perpendicular to B
        for xn = 1:xR        
%         for xn = 5
%         for xn = Rnxn(2)
            x = xn*dx;

            GRx = x/stheta + b*(-R*stheta + x*ctheta)^2;
%             GRx = 2*x*ctheta+R*(stheta^2-ctheta^2)+b*stheta*(-R*ctheta+x*stheta)^2;
            
            if GRx > 0

                vb = 0;
                deltavb = 0;
                Nb = 0;
                integralF = 0;
                fi = 0;

                %%  new numerical method with higher resolution around vc                              
                Xp = Xt + x*stheta + R*ctheta;
                Yp = Yt + x*ctheta - R*stheta;
                               
%                 if R >= 0 
%                     vc = 0.5*vd*(R/x + ctheta/stheta + sqrt((R/x - ctheta/stheta)^2 + 1/(b*x*stheta^3)));
%                 else
%                     vc = 0.5*vd*(R/x + ctheta/stheta - sqrt((R/x - ctheta/stheta)^2 + 1/(b*x*stheta^3)));
%                 end
                if R>=0
                    Ys = Yp - sqrt(Yp^2 + (Xp-as)/b);
                else
                    Ys = Yp + sqrt(Yp^2 + (Xp-as)/b);
                end
                Xs = as - b*Ys^2;
                
                vc = sign(R)*vsw*abs(1/(ctheta - 2*b*Ys*stheta) - ctheta);                
                
                if abs(vc) < vbmin || abs(vc) > vmax
                    vc = 0;
                    vbRx(Rn,xn) = 0;
                else                   
%                     vcRx(Rn,xn) = vc;
                    stop = 0;
                    Yb = Ys;
%                     Rb = Rc;
                    vel_b = vc;
                    vel_a = vc;                    

                    vel_a0 = 0;
                    vel_b0 = 0;
                    
                    iter = sign(R)*iter0;
                    
                    if not(B2B1step)
                        Y1 = Ys;
                        thetbn1 = acosd((ctheta - 2*b*Y1*stheta)./sqrt(1+(2*b*Y1).^2+(2*b*Z)^2));
                        c2thetabnmin1 = cosd(thetbn1).^2;
                        s2thetabnmin1 = sind(thetbn1).^2;

                        % Kivelson & Russel
                        gam = 5/3;
                        thetvn1 = acosd(1./sqrt(1+(2*b*Y1).^2+(2*b*Z)^2));

                        Mamin1 = Matot/sqrt(1+(2*b*Y1).^2+(2*b*Z)^2);
                        M_s1 = vsw./sqrt(1+(2*b*Y1).^2+(2*b*Z)^2)/cs;

                        % Alfven speed and fast mode mach number
    %                             B0 = 6.4e-9;
    %                             va=sqrt(B0^2/mu0/mi/Ne);
                        va = vsw/Matot;
    %                                     vfast = sqrt(0.5*((va^2+cs^2 + sqrt((va^2 + cs^2).^2-4*va^2*cs^2.*ctheta^2))));
    %                                     M_f = vsw./sqrt(1+(2*b*Y).^2+(2*b*Z)^2)/vfast;

    %                                     if sum(M_f < 1) > 0
    %                                         M_f_constraint = M_f_constraint + 1;
    %                                         B2B1 = ones(1,step+1);
                        if Mamin1 > 1

                            Ma1 = Mamin1;
                            Ms1 = M_s1;

                            na1 = -((gam-1)/(gam+1) + (gam+(gam+2)*cosd(thetbn1).^2)./((gam+1)*Ma1.^2) + 2./((gam+1)*Ms1.^2));
                            nb1 = (gam*(1+cosd(thetbn1).^2) - 2 + cosd(thetbn1).^2.*((gam+1)./Ma1.^2 + 4./Ms1.^2))./((gam+1)*Ma1.^2);
                            nc1 = -cosd(thetbn1).^2./((gam+1)*Ma1.^4).*(gam - 1 + 2*cosd(thetbn1).^2./Ms1.^2);

                            np1 = (3*nb1 - na1.^2)/3;
                            nq1 = (2*na1.^3)/27 - na1.*nb1/3 + nc1;
                            nD1 = (np1/3).^3 + (nq1/2).^2;                                        

                            nu1 = (-nq1/2 + sqrt(nD1)).^(1/3);
                            nv1 = (-nq1/2 - sqrt(nD1)).^(1/3);                            

                            % density jump across shock
                            n2n11 = 1./(nu1 + nv1 - na1/3);                                        
                        else
                            % density jump across shock
%                             n2n11 = 1./((gam-1)/(gam+1)+2/(gam+1)./M_s1.^2);
                            Mamin1 = M_s1;
%                             n2n11 = ((gam+1)*M_s1^2)/((gam-1)*M_s1^2+2);
                            n2n11 = 6/(1+5*M_s1^(-2));
                        end
                        % tangential B-field jump
                        Bt2Bt1min1 = n2n11.*(Mamin1.^2-c2thetabnmin1)./(Mamin1.^2-n2n11.*c2thetabnmin1);
                        B2B1 = sqrt(c2thetabnmin1 + s2thetabnmin1.*Bt2Bt1min1.^2);
                    end

                    while stop == 0

                        Yb = Yb + iter;
%                         if Yb ~= Ys                        

                            Xb = as - b*Yb.^2;
                            Rb = -(Yb-Yt)*stheta+(Xb-Xt)*ctheta;
                            xb = (Yb-Yt)*ctheta+(Xb-Xt)*stheta;

                            vel_b = vd*(R-Rb)/(x-xb);

                            if abs(vel_b) > c 
                                if abs(iter) > itermin 
                                    Yb = Yb - iter;
                                    iter = iter/10;
                                else
                                    vel_a = vc;
                                    vel_b = vc;
                                    stop = 1;
                                end  
                            elseif abs(vel_b) > abs(vc) %&& abs(vel_b) > abs(vel_b0)               
                                step = 1000;

                                if R >= 0
                                    h = (Yb - Ys)/step;
                                    Y = Ys:h:Yb;
    %                                 h = (Ys - Yb)/step;
    %                                 Y = Yb:h:Ys;

                                else
    %                                 h = (Yb - Ys)/step;
    %                                 Y = Ys:h:Yb;
                                    h = (Ys - Yb)/step;
                                    Y = Yb:h:Ys;
                                end

                                X = as - b*Y.^2;

                                Rs = -(Y-Yt)*stheta+(X-Xt)*ctheta;
                                xs = (Y-Yt)*ctheta+(X-Xt)*stheta;
                                v = vd*(R-Rs)./(x-xs);

%                                 thetbn = acosd((ctheta - 2*b*Y*stheta)./sqrt(1+(2*b*Y).^2+(2*b*Z)^2));

                                vpht = vsw*(2*b*Y*stheta*ctheta + stheta^2)./(ctheta - 2*b*Y*stheta);                                                
    %                             vpht = vd*tand(thetbn);

                                vparai2 = (2.*vpht-v).^2;
    %                             vparai2 = (vpht-v).^2;

%                                 if B2B1step
%                                     thetbn = acosd((ctheta - 2*b*Y*stheta)./sqrt(1+(2*b*Y).^2+(2*b*Z)^2));
%                                     c2thetabnmin = cosd(thetbn).^2;
%                                     s2thetabnmin = sind(thetbn).^2;
% 
%                                     % Kivelson & Russel
%                                     gam = 5/3;
%                                     thetvn = acosd(1./sqrt(1+(2*b*Y).^2+(2*b*Z)^2));
% 
%                                     Mamin = Matot./sqrt(1+(2*b*Y).^2+(2*b*Z)^2);
%                                     M_s = vsw./sqrt(1+(2*b*Y).^2+(2*b*Z)^2)/cs;
% 
%                                     % Alfven speed and fast mode mach number
%         %                             B0 = 6.4e-9;
%         %                             va=sqrt(B0^2/mu0/mi/Ne);
%                                     va = vsw/Matot;
% %                                     vfast = sqrt(0.5*((va^2+cs^2 + sqrt((va^2 + cs^2).^2-4*va^2*cs^2.*ctheta^2))));
% %                                     M_f = vsw./sqrt(1+(2*b*Y).^2+(2*b*Z)^2)/vfast;
% 
% %                                     if sum(M_f < 1) > 0
% %                                         M_f_constraint = M_f_constraint + 1;
% %                                         B2B1 = ones(1,step+1);
%                                     if sum(Mamin < 1) == 0
% 
%                                         Ma = Mamin;
%                                         Ms = M_s;
% 
%                                         na = -((gam-1)/(gam+1) + (gam+(gam+2)*cosd(thetbn).^2)./((gam+1)*Ma.^2) + 2./((gam+1)*Ms.^2));
%                                         nb = (gam*(1+cosd(thetbn).^2) - 2 + cosd(thetbn).^2.*((gam+1)./Ma.^2 + 4./Ms.^2))./((gam+1)*Ma.^2);
%                                         nc = -cosd(thetbn).^2./((gam+1)*Ma.^4).*(gam - 1 + 2*cosd(thetbn).^2./Ms.^2);
% 
%                                         np = (3*nb - na.^2)/3;
%                                         nq = (2*na.^3)/27 - na.*nb/3 + nc;
%                                         nD = (np/3).^3 + (nq/2).^2;                                        
% 
%                                         nu = (-nq/2 + sqrt(nD)).^(1/3);
%                                         nv = (-nq/2 - sqrt(nD)).^(1/3);                            
% 
%                                         % density jump across shock
%                                         n2n1 = 1./(nu + nv - na/3);
% 
%                                         % tangential B-field jump
%                                         Bt2Bt1min = n2n1.*(Mamin.^2-c2thetabnmin)./(Mamin.^2-n2n1.*c2thetabnmin);
%                                         B2B1 = sqrt(c2thetabnmin + s2thetabnmin.*Bt2Bt1min.^2);
%                                     else
%                                         % density jump across shock
%                                         n2n1 = 1./((gam-1)/(gam+1)+2/(gam+1)./M_s.^2);
% 
%                                         % tangential B-field jump
%                                         Bt2Bt1min = n2n1.*(Mamin.^2-c2thetabnmin)./(Mamin.^2-n2n1.*c2thetabnmin);
%                                         B2B1 = sqrt(c2thetabnmin + s2thetabnmin.*Bt2Bt1min.^2);
%                                     end                                                                                                       
%                                 end

                                if sum(B2B1 < 1) > 0
                                    B2B1 = ones(1,step+1);
                                end
        %                         B2B1Rx(Rn,xn) = B2B1(1);

                                deltaphi = 2/me*3*kbolt*Te*(B2B1-1);

                                vlc2i = ((v-vpht).^2+deltaphi)./(B2B1-1);

                                gamma_F = nhalo*gamma(kappa+1)/gamma(kappa-0.5)/(ve*sqrt(2*pi*kappa^3));
                                gkf = gamma_F*(2*kappa*ve^2)^kappa;


                                Fi = gkf*(2.*kappa.*ve.^2+v.^2).^(-kappa);                                  
                                
                                Fr = gkf*(2*kappa*ve^2 + vlc2i+vparai2).^(-kappa);

                                F = Fi + Fr;                                               
        %                         
        %                         FZ(Zn,:) = F;
                                if F(step/2+1+sign(R)*step/2) ~= max(F)
                                    if R >= 0
                                        vel_a = sign(vel_b)*sqrt((gkf/F(step+1))^(1/3)-2*kappa*ve^2);
                                    else
                                        vel_a = sign(vel_b)*sqrt((gkf/F(1))^(1/3)-2*kappa*ve^2);
                                    end

                                    if abs(vel_a) < abs(vc)
                                        cva = 2*kappa*ve^2;

                                        fivel_a = vel_a./(4*cva*(vel_a.^2+cva).^2) + 3*vel_a./(8*cva^2*(vel_a.^2+cva)) + 3*atan(vel_a./sqrt(cva))./(8*cva^(2.5));
                                        fivc = vc./(4*cva*(vc.^2+cva).^2) + 3*vc./(8*cva^2*(vc.^2+cva)) + 3*atan(vc./sqrt(cva))./(8*cva^(2.5));

                                        fi = gkf*abs(fivc-fivel_a);

                                        Fi_a = gkf*(2.*kappa.*ve.^2+vel_a.^2).^(-kappa);

                                        A = abs(vel_a-vc).*Fi_a-fi;

    %                                     simpsonF = (F(1)+F(step+1)+4*sum(F(2:2:step))+2*sum(F(3:2:step-1)))*h/3;
    %                                     simpsonerrorF = (F(1)+F(step+1)+4*sum(F(3:4:step-1))+2*sum(F(5:4:step-3)))*2*h/3;
    %                                     integralF = simpsonF + (simpsonF-simpsonerrorF)/15;
                                        integralF = F(1)*abs(v(1)-v(2))/2;
                                        integralF = integralF + F(step+1)*abs(v(step+1)-v(step))/2;
                                        integralF = integralF + sum(F(2:step).*abs(v(3:step+1)-v(1:step-1))/2);
%                                         for i = 2:step                                        
%                                             integralF = integralF + F(i)*abs(v(i+1)-v(i-1))/2;
%                                         end

                                        if R >= 0
                                            B = integralF - abs(vel_b-vc).*F(step+1);                                       
                                        else
                                            B = integralF - abs(vel_b-vc).*F(1);                        
                                        end

                                        mini = B - A;

                                        if mini > 0
                                            if abs(iter) > itermin
                                                Yb = Yb - iter;
                                                iter = iter/10;
                                            else
                                                stop = 1;
                                            end
                                        else                                    
                                            if vel_a == vel_a0
                                                stop = 1;
                                            else
                                                vel_a0 = vel_a;
                                                vel_b0 = vel_b;
                                            end
                                        end
                                   else
                                        if abs(iter) > itermin
                                            Yb = Yb - iter;
                                            iter = iter/10;                                
                                        else
                                            vel_a = vc;
                                            vel_b = vc;
                                            stop = 1;
                                        end
                                    end
                                else
                                    vel_a0 = vel_a;
                                    vel_b0 = vel_b;
                                end
                            else
                                if abs(iter) > itermin
                                    Yb = Yb - iter;
                                    iter = iter/10;
                                else
                                    vel_a = vc;
                                    vel_b = vc;
                                    stop = 1;
                                end
                            end
%                         end
                    end                                        
                                       
                    vb = 0.5*(vel_a+vel_b);
                    deltavb = 0.5*abs(vel_b - vel_a);
                    Nb = (integralF + fi);
%                     Nb2 = gkf*(2.*kappa.*ve.^2+vel_a^2).^(-kappa)*abs(vel_b - vel_a);
                    
%                     vcRx(Rn,xn,Zn) = vc;
                    vbRx(Rn,xn,Zn) = vb;
%                     [arb,indF] = max(F);
%                     vbRx2(Rn,xn,Zn) = v(indF);
                    NbRx(Rn,xn,Zn) = Nb;
                    dvbRx(Rn,xn,Zn) = deltavb;
%                     dvbRx2(Rn,xn,Zn) = vel_b-vc;                    
                    Wfr(Rn,xn,Zn) = Nb*me*abs(vb)*deltavb/3;    
                    
                    %for comparison with line of flight
                    l(Rn,xn,Zn) = sqrt((Xp-Xs)^2+(Yp-Ys)^2);
                                        
%                     % Schmidt comparison
%                     schmidt1 = 2*kappa*ve^2;
% %                     schmidt1 = ve^2;
%                     schmidt2 = schmidt1*(B2B1-1);
% %                     schmidt3 = 1+vc^2/schmidt1+schmidt2;
%                     schmidt4 = 1+vc^2/schmidt1+deltaphi/schmidt2;
%                     schmidt5 = sqrt(1+((1+vc^2/schmidt1)/schmidt4)^(kappa+1))-1;
%                     schmidt6 = (schmidt4/(1+vc^2/schmidt1))^(kappa+1);
%                     a_v = vc - ve^2/vc*schmidt4*schmidt5;
%                     b_v = vc + ve^2/vc*schmidt4/(1+schmidt6)*(1-schmidt6*schmidt5);
%                     
%                     %new calc
%                     %cva = 2*kappa*ve^2
%                     a_v2 = vc - cva/vc*(1+vc^2/cva+deltaphi/(cva*(B2B1-1)))*(sqrt(1+((1+vc^2/cva)/(1+vc^2/cva+deltaphi/(cva*(B2B1-1))))^(kappa+1))-1);
%                     b_v2 = vc + cva/vc*(1+vc^2/cva+deltaphi/(cva*(B2B1-1)))/(1+((1+vc^2/cva+deltaphi/(cva*(B2B1-1)))/(1+vc^2/cva))^(kappa+1))*(1-((1+vc^2/cva+deltaphi/(cva*(B2B1-1)))/(1+vc^2/cva))^(kappa+1)*(sqrt(1+((1+vc^2/cva)/(1+vc^2/cva+deltaphi/(cva*(B2B1-1))))^(kappa+1))-1));
%                     
%                     vb2 = 0.5*(a_v+b_v);
%                     deltavb2 = 0.5*abs(b_v - a_v);
%                     Nb3 = gkf*(2.*kappa.*ve.^2+a_v^2).^(-kappa)*abs(b_v - a_v);
% %                     
% %                     vcRx(Rn,xn,Zn) = vc;
%                     vbRx(Rn,xn,Zn) = vb2;
% %                     [arb,indF] = max(F);
% %                     vbRx2(Rn,xn,Zn) = v(indF);
%                     NbRx(Rn,xn,Zn) = Nb3;
%                     dvbRx(Rn,xn,Zn) = deltavb2;
% %                     dvbRx2(Rn,xn,Zn) = vel_b-vc;                    
%                     Wfr(Rn,xn,Zn) = Nb3*me*abs(vb2)*deltavb2/3;    

                    %                     figure()
%                     hold on
%                     semilogy(v,F)                                                            
%                     plot(v,F)
                end        
            end   
        end
%         toc
    end
    toc
%     if M_f_constraint > 0
    if Z > 100*Rp
        loop = false;
    end
end
    
    %% Efficiencies
    alphas = 1;
    beta = 1/3;
    gammaLS = temp4gammaLS./(alphas*vbRx.^2);
    zetaF = exp(-temp4zetaF.*(abs(vbRx)./(beta*dvbRx)).^2.*((vcrit/ve-abs(vbRx)./ve).^2));
    temp4zetaH = ve*beta*dvbRx*sqrt(2)./vbRx.^2;
    zetaH = sqrt(pi/6)*beta/2*(erf(real(erftemp1./temp4zetaH))+erf(real(erftemp2./temp4zetaH)));
%     phiFt = temp4phiF.*zetaF.*gammaLS;
    phiFt = temp4phiF.*zetaF.*gammaLS/2*3.*abs(vbRx)./dvbRx;
    
    phiH = temp4phiH.*zetaH.*abs(vbRx)/c;
    
    testdiff3D

%     %% power density
% 
%     % Differentiation via backward differencing (ie from upstream regions)
%     dWvbdR = zeros(size(vbRx));
%     dWdx = zeros(size(vbRx));
% for Zn = 1:size(vbRx,3)
%     for Rn = 1:RR/2-2
%         dWvbdR(Rn,:,Zn)= -(Wfr(Rn,:,Zn).*vbRx(Rn,:,Zn)- Wfr(Rn+2,:,Zn).*vbRx(Rn+2,:,Zn))/2/dR;
%     end
%     for Rn = RR/2+3:RR
%         dWvbdR(Rn,:,Zn)= (Wfr(Rn,:,Zn).*vbRx(Rn,:,Zn)- Wfr(Rn-2,:,Zn).*vbRx(Rn-2,:,Zn))/2/dR;
%     end
%     dWvbdR(Rn/2-1,:,Zn)= -(Wfr(RR/2-1,:,Zn).*vbRx(RR/2-1,:,Zn) - Wfr(RR/2,:,Zn).*vbRx(RR/2,:,Zn))/dR;
%     dWvbdR(RR/2,:,Zn)= dWvbdR(RR/2-2,:,Zn)+dWvbdR(RR/2-3,:,Zn)-dWvbdR(RR/2-1,:,Zn);
%     dWvbdR(RR/2+1,:,Zn)= dWvbdR(RR/2+3,:,Zn)+dWvbdR(RR/2+4,:,Zn)-dWvbdR(RR/2+2,:,Zn);
%     dWvbdR(RR/2+2,:,Zn)= (Wfr(RR/2+2,:,Zn).*vbRx(RR/2+2,:,Zn) - Wfr(RR/2+1,:,Zn).*vbRx(RR/2+1,:,Zn))/dR;
%     for xn=3:xR
%         dWdx(:,xn,Zn)= (Wfr(:,xn,Zn)-Wfr(:,xn-2,Zn))/2/dx;
%     end
% end
%     dWdx(:,2,:) = (Wfr(:,2,:)-Wfr(:,1,:))/dx;
%     dWdx(:,1,:) = dWdx(:,3,:)+dWdx(:,4,:)-dWdx(:,2,:);
% 
%     Wflux = dWvbdR + vd*dWdx;
% % 
% %     inc = -1;
% % for Rn = 1:RR
% %     if Rn>RR/2
% %         inc = 1;
% %     end
% %     dRn = inc+40*ceil(-1+2*(Rn-1)/(RR-1));
% %     Rleft = min(Rn,Rn-dRn);
% %     if Rleft < 1
% %         Rleft = 1;
% %     end        
% %     Rright = max(Rn,Rn-dRn);
% %     for xn = 1:xR
% %         if abs(vbRx(Rn,xn))>0 && abs(vbRx(Rn,xn))<vbmin
% %             break
% %         end
% %         if Wflux(Rn,xn) ~= 0
% %             dxn = 5+2*ceil(100*xn/xR);
% %             xleft = xn-dxn;
% %             if xleft < 1
% %                 xleft = 1;
% %             end
% %             Wflux_smooth(Rn,xn)=sum(sum(Wflux(Rleft:Rright,xleft:xn)))/(size(Wflux(Rleft:Rright,xleft:xn),1)*size(Wflux(Rleft:Rright,xleft:xn),2));
% %         end
% %     end
% % end
% %     
% % Wflux = Wflux_smooth;
% % phiFt(isnan(phiFt)) = 0;
% %% emissivities
% % power from Kuncic paper and LMC
% jF = phiFt.*Wflux/coneF;
% jF(isnan(jF)) = 0;
% jH = phiH.*Wflux/coneH;
% jH(isnan(jH)) = 0;
% % jF_nlt2 = phiFt2.*Wflux/coneF;
% % LMC efficiency
% etae = 0.3*(3*ve/c)^4/(sqrt(0.0167)*pi)*exp(-Uc^2)/(sqrt(pi)*Uc);
% jF_lmc = etae.*Wflux/coneF.*(abs(phiFt)>0); 
% jF_lmc(isnan(jF_lmc)) = 0;
%     