ag/cross section fix #428
ag/cross section fix #428
Conversation
| @@ -333,19 +340,15 @@ def cross_section(self, phi, thetas=None): | |||
| single curve. | |||
There was a problem hiding this comment.
Can you add a description of the arguments, i.e.
inputs
-------
phi: float, cylindrical angle between [0, 2pi] defining the cross section.
thetas: optional, Union[int, arrayLike, None]. Poloidal angle values at which to evaluate the surface on the cross section. If thetas is an integer, then the same number of evenly spaced points are used. Alternatively, thetas can be an array of poloidal angles values between [0,1].
| @@ -333,19 +340,15 @@ def cross_section(self, phi, thetas=None): | |||
| single curve. | |||
| """ | |||
|
|
|||
| # phi is assumed to be between [-pi, pi], so if it does not lie on that interval | |||
| # phi is assumed to be between [0, 2*pi], so if it does not lie on that interval | |||
There was a problem hiding this comment.
Can we make phi and theta follow the same convention, i.e. both in [0,1]?
| @@ -333,19 +340,15 @@ def cross_section(self, phi, thetas=None): | |||
| single curve. | |||
| """ | |||
|
|
|||
| # phi is assumed to be between [-pi, pi], so if it does not lie on that interval | |||
| # phi is assumed to be between [0, 2*pi], so if it does not lie on that interval | |||
| # we shift it by multiples of 2pi until it does | |||
There was a problem hiding this comment.
Should we be correcting the inputs or expecting phi to be within a range? If we expect phi in [0, 1] or [0, 2pi] then we users should coerce their inputs to [0, 1] or [0, 2pi]? We would then can replace this code with assert (0 <= phi) & (phi <= 1), "phi must be in [0,1]. If we make this change, we should look to make sure we aren't breaking usage in other functions.
Alternatively, in the docstring we could specify that any value of phi is acceptable. Then this code would stay.
What do you think?
| cyl_phi_right = cyl_phi[1, :] | ||
|
|
||
| phi = phi*np.ones(varphi_left.size) | ||
| def put_angle_between_left_and_right(angle, left_bound, right_bound): |
There was a problem hiding this comment.
can you add some description. function name says angle, is that phi or varphi?
| phi = phi - 2. * np.pi | ||
| if phi < -np.pi: | ||
| phi = phi + 2. * np.pi | ||
|
|
There was a problem hiding this comment.
if phi is the cylindrical phi for our surface, can we just return the points from gamma?
| angle = np.arctan2(gamma[:, 1], gamma[:, 0]) | ||
| angle = np.where(angle < 0, angle+2*np.pi, angle) | ||
| shifted_angle = put_angle_between_left_and_right(angle, left_bound, right_bound) | ||
| return shifted_angle | ||
|
|
||
| def bisection(phia, a, phic, c): |
There was a problem hiding this comment.
can we use scipy bisection or root finding?
| cyl_phi_right = cyl_phi[idx_right, np.arange(idx_right.size)] | ||
|
|
||
| # this function converts varphi to cylindrical phi, ensuring that the returned angle | ||
| # lies between left_bound and right_bound. | ||
| def varphi2phi(varphi_in, left_bound, right_bound): |
There was a problem hiding this comment.
describe a little here as well
| k = np.ceil((left_bound-angle)/(2*np.pi)) | ||
| shifted_angle = angle + 2*np.pi * k | ||
| if not np.all((left_bound <= shifted_angle) & (shifted_angle <= right_bound)): | ||
| import ipdb;ipdb.set_trace() |
There was a problem hiding this comment.
remove debugging code?
| cyl_phi = np.where(cyl_phi < 0, cyl_phi + 2*np.pi, cyl_phi) | ||
| cyl_phi[1, :] += 2*np.pi # second row is the same as the first row, but shifted by 2pi | ||
|
|
||
| varphi_left = varphigrid[0, :] |
There was a problem hiding this comment.
Are we trying to determine bisection bounds here?
| @@ -425,7 +411,8 @@ def bisection(phia, a, phic, c): | |||
| err = np.max(np.abs(a - c)) | |||
| b = (a + c) / 2. | |||
| return b | |||
|
|
|||
|
|
|||
| phi = put_angle_between_left_and_right(phi, cyl_phi_left, cyl_phi_right) | |||
| # bisect cyl_phi to compute the cross section | |||
| sol = bisection(cyl_phi_left, varphi_left, cyl_phi_right, varphi_right) | |||
There was a problem hiding this comment.
just trying to understand here: it looks like, for each theta, we are doing a bisection over varphi, so that varphi2phi(varphi) = phi. Why do we need to know cyl_phi_right and cyl_phi_left?
I have used the
Surface.cross_sectionalgorithm on the entire QUASR database, and notice that very infrequently, the algorithm fails, see belowI rewrote the algorithm, making it simpler, to fix the bug. Rerunning the new algorithm on the entire dataset, I can't find any incorrect cross sections anymore.
Note that the jagged cross section where
varphi=0is because it uniformly samples the cross section in the Boozer theta angle, which may be insufficient as you can see.