definition of beam radius in gaussian intensity
BrianBostwick opened this issue · 6 comments
I think that these are not properly defined,
Line 71 in 79147a4
Line 72 in 79147a4
and
Line 105 in 79147a4
Line 106 in 79147a4
should they not be
let std = e_radius / 2.0_f64.powf(0.5);
let power = std::f64::consts::PI * std.powi(2) * peak_intensity;
When they are changed to this above however the unit tests failures are:
integration_tests::rate_equation::tests::single_beam_scattering_rates_v_detuning
integration_tests::rate_equation::tests::single_beam_scattering_rates_v_intensity
If someone would take a look and see if they agree with this change or please explain why this is correct?
For some clarity I'll redefine some relations here.
Taking an infinite Rayleigh range, the gaussian intensity in the program is defined as:
power / PI / beam.e_radius.powf(2.0)
* EXP.powf(
-distance_squared
/ beam.e_radius.powf(2.0),
)
as (the 1/rayleigh_range terms go to zero for infinite range).
Rewriting this in symbolic terms:
intensity, I = P / (pi r_0^2) exp (-r^2/r_0^2)
Re-deriving the above
We define the gaussian intensity as:
I = I_0 exp (-r^2 / r_0^2)
where the constant r_0 is defined as the 1/e-radius of the intensity (NOT the 1/e^2 radius, as is used for a beam waist definition - see documentation for comment):
I = I_0 exp (-1) for r = r_0
We relate the factor I_0 to the power of the beam using the 2D gaussian integral:
P = \int I dx^2
P = I_0 \int_0^{2pi} d \theta \int_0^\infty r exp(-r^2 / r_0^2) dr
P = I_0 [2 \pi] [r_0^2 / 2]
P = \pi r_0^2 I_0
Rearranging for I_0, we recover the power / PI / beam.e_radius.powf(2.0) factor in the code.
So the code looks good to me, but I might have made a mistake. Can I ask you to double check you are not confusing the 1/e-radius with the beam waist? The beam waist \omega_0 is defined as the diameter for which I = I_0 exp(-2), which can be rearranged to find:
I = I_0 exp(-(\omega/2)^2/r_0^2) = exp(-2)
\rightarrow 2 r_0^2 = (\omega_0 / 2)^2
\omega_0 = 2 \sqrt{2} r_0
(Also I feel slightly mad, maybe my brain has already run out today. In your post above, is there anything different between your lines of what they are previously and what they should be? I can't see the difference!)
Yes, I agree. For some reason when I look at this it makes me think that std is the e^2 radius (which is why I suggested the change / -> *) and the 2 pre-factor should be a 1/2 then. I think I was confused by the factors of 2 in there. Is there a reason why they have been added? Since,
power = (e_radius)^2 * pi * I
Why not just put that in there directly.
The choice of constant is somewhat arbitrary (e.g. whether parameterized as e_radius or beam_waist). However, I would be very hesitant to change the definition of e_radius now, because it would alter the behavior without warning (e.g. in programs where people have set e_radius themselves - such as the unit tests, so it's good those broke because they represent that some part of the program has changed for equivalent inputs).
The best approach would be to improve the API, e.g. by adding a constructor for GaussianBeam that could allow you to define the beam waist, e.g. GaussianBeam::create().with_power(x).with_beam_waist(w0). This would loop back on #41 .
If you really want to change the parameterization, I think e_radius should be renamed to something clearer (because the variable's intent has changed, but also to force a breaking change so that people are aware the meaning of the quantity has changed). However, my preference remains to not change it and instead improve the public facing API.
Apologies if misunderstanding you!
No that makes sense, I agree, I don't think it's worth going through and changing everything now. What did you have in mind when it comes to
My comment had more to do with why there is a value called "std" it seems like it is just adding a factor of sqrt two to remove the two in the next line. It was confusing because it makes it look like e_radius is actual e^2_radius when its not. So I wanted to know if "std" served a prepose other then the removing the twos and what it was? If it does not I was going to propose a change to following and removing "std"
` let power = std::f64::consts::PI * e_radius.powi(2) * peak_intens
I don's know what factory structs are but I can take a look at what you are thinking and see if I can implement those changes to GaussianBeam. Should we continue that conversation in #41 then? and Ill close this issue.
My comment had more to do with why there is a value called "std" it seems like it is just adding a factor of sqrt two to remove the two in the next line. It was confusing because it makes it look like e_radius is actual e^2_radius when its not. So I wanted to know if "std" served a prepose other then the removing the twos and what it was? If it does not I was going to propose a change to following and removing "std"
It does look like there is some redundancy here and you are right that it could be simplified. I guess it's defining the standard deviation and then using a textbook result for the area under a gaussian curve.
Let's close the issue and discuss API changes in the other issue.