rsun=6.96×10^10
Rstar = 19 rsun
Rmin = 1.5*Rstar
uo = Rstar/Rmin
vinf = 2.25×10^8
beta = 1.0
mdotmsunyr = 8.×10^(-6)
mdot = 1.99×10^(33)* mdotmsunyr/(3.15×10^7)
c = 1.
q = 0.

v[r_]:=vinf*(1.-Rstar/r)^beta
Plot[v[r], {r,Rstar,10Rstar}]

rho[r_] := (mdot/(4 π r^2  v[r]))

Plot[rho[r], {r,Rstar, 10 Rstar}]

w[r_] := v[r]/vinf
Plot[w[r],{r,Rstar,10 Rstar}];

u[r_] := Rstar/r
v[u_] := vinf (1.-u)^beta

Plot[v[u],{u,1,0}];

w[u_] := v[u]/vinf

Plot[w[u],{u,1,0}];

f[r_] := (r/Rmin)^(-q)

Plot[f[r],{r,Rstar, 10 Rstar}];

f[u_] := (u*(Rmin/Rstar))^q

Plot[f[u], {u,1,0}, AxesLabel->Automatic];

emissivity[x_,r_,μ_] := c rho[r]^2 DiracDelta[x+(μ*w[r])]

c1 = (c*mdot^2/2/vinf^2/Rstar)

1.9076871167364177`*^12

** OK, ALL THE VARIABLES ARE SET, AND FUNCTIONS DEFINED -- EXCEPT FOR OPTICAL DEPTH **

** NOW CALCULATE PROFILE(S) FOR A COMPLETELY OPTICALLY THIN LINE **
** FOR NOW WITHOUT OCCULTATION **

usubx[x_,beta_] := Min[uo,1.-Abs[x]^(1/beta)]

xlum[x_,beta_] := c1*(∫_0^(usubx[x,beta]))((f[u])/((1-u)^(3*beta)))\[DifferentialD]u

Plot[{xlum[x,1.],xlum[x,3.]},{x,-1,1}]

** NOTE TWO TEST VALUES OF BETA ABOVE, BUT BELOW HERE, I NORMALIZE ONE PROFILE.  I NEED A WAY TO NORMALIZE EACH PROFILE (FOR DIFFERENT VALUES OF BETA) RELATIVELY AUTOMATICALLY **

** FOR NOW I JUST NORMALIZE ONE VALUE (WHICH APPEARS TO BE THE DEFAULT VALUE, SET AT THE TOP OF THIS NOTEBOOK TO ZERO) **

maxxlum =First[ FindMaximum[{xlum[x,beta],-1≤x≤1},{x,-1}]]

normxlum[x_,beta_] := xlum[x,beta]/maxxlum

Plot[normxlum[x,beta],{x,-1,1}, AxesLabel->Automatic]

** OCCULTATION **

** OPTICAL DEPTH **