function [t,L]=L_of_M_vmax_kappa_xs_ts(M,v,kappa,xs,ts,Q0,tmax,Nt,nmax)
%returns Luminosity L (erg/s) and time t (seconds) for an ejecta with mass M (grams)
%with uniform density, with the outer edge expanding with velocity v
%(cm/s), with a deposition that is uniform within inner fraction xs of the
%radius, and total deposition as function of time given by  Q=Q0*exp(-t/ts), where ts is the
%decay time in seconds and Q0 is the deposition at t=0 in erg/cm^3. tmax is the time to which to calculate (seconds)
%Nt is the number of time steps (uniform in t^2). nmax is number of terms
%to include in the expansion
% example use:
%M=2e33; v=1e9; kappa=0.1; ts=10*86400; Q0=5e43; xs=0.3;
%tmax=100*86400; Nt=10000; nmax=100; 
%[t,L]=L_of_M_vmax_kappa_xs_ts(M,v,kappa,xs,ts,Q0,tmax,Nt,nmax);



c=2.99792458e10; %cm/s
rhot=M/(4*pi*v^3/3); %g s^3 cm^-3
D=c/(3*kappa*rhot); %cm^2 s^-4 
t_unit_tilde=v^2/D;  %s^2
E_unit_tilde=Q0*ts^2; %erg*s (=\int Qtdt)
L_unit_tilde=E_unit_tilde/t_unit_tilde; %erg/s

tmax_tilde=0.5*tmax^2; % s^2
tmax_tilde_norm=tmax_tilde/t_unit_tilde; %dimensionless
 
ts_tilde=0.5*ts^2; % s^2
ts_tilde_norm=ts_tilde/t_unit_tilde; %dimensionless
[t_tilde_norm,L_tilde_norm]=L_of_t_xs_ts_nmax(xs,ts_tilde_norm,nmax,tmax_tilde_norm,Nt);
t_tilde=t_tilde_norm*t_unit_tilde; %s^2
t=sqrt(2*t_tilde); %s
L_tilde=L_tilde_norm*L_unit_tilde; %erg/s
L=L_tilde; %erg/s
end


function [t,L,Ldelta]=L_of_t_xs_ts_nmax(xs,ts,nmax,tmax,Nt)
    %Luminosity for diffusion in a uniform static ball of radius R=1 and diffusion coefficient
    %D=1, with energy deposited uniformly over an inner ball with radius xs and decaying 
    %as a function of time as  Q(t)=exp(-sqrt(t/ts))/ts/2 (total integral equal unity) 
    %nmax is the number of terms to include in the expansion. 
    t=linspace(0,tmax,Nt);
    dt=tmax/(Nt-1);
    Ldelta=L_of_t_xs_nmax(t,xs,nmax);
    L=t*0;  Nt=length(t);
    tic
    if ts/100>dt
        for lt=1:Nt
            tt=(lt-0.5)*dt;
            L(lt:end-1)=L(lt:end-1)+0.5*(Ldelta(1:end-lt)+Ldelta(2:end-lt+1))*exp(-sqrt(tt/ts))*dt/ts/2;
        end
    else  
        dtint=ts/40;
        tint=dtint/2 : dtint : 100*ts ;
        Ntint=length(tint);
        for lt=1:Ntint
            tnow=tint(lt);
            L=L+interp1(t,Ldelta,t-tnow,'linear',0)*exp(-sqrt(tnow/ts))*dtint/ts/2;
        end
    end
    t=t(1:end-1); L=L(1:end-1);
end

function [L,fout]=L_of_t_xs_nmax(t,xs,nmax)
    %Luminosity and accumulated luminosity for diffusion in a uniform static ball of radius R=1 and
    %diffusion coefficient D=1, with an impulse deposition of unit energy spread
    %uniformly over an inner ball of radius xs. nmax is the number of terms
    %to include in the expansion. 
    L=t*0; fout=t*0+1;
    n=[1:nmax];
    if xs==0
        An=2*n.^2*pi^2.*(-1).^(n+1);  
    else
        An=6/xs^3*(-1).^(n+1).*(sin(pi*n*xs)./n/pi-xs*cos(n*pi*xs));
    end

    for n=1:nmax
        Ln=An(n)*exp(-n^2*pi^2*t);
        L=L+Ln;
        fout=fout-Ln/(n^2*pi^2); % int_0^t L(t')dt'
    end

    ind0=find(t==0);
    L(ind0)=0; fout(ind0)=0;
end