improve documentation for large function

[hematthi]

EchelleCCFs.jl/src/ccf/calc_ccf.jl

Lines 304 to 413 in a639217

	""" `project_mask!( output, λs, ccf_plan; shift_factor )`
	Compute the projection of the mask onto the 1D array of wavelengths (λs) at a given shift factor (default: 1).
	The mask is computed from the ccf_plan, including a line_list and mask_shape (default: tophat).
	Assumes plan.line_list is already sorted.
	"""
	function project_mask!(projection::A2, λs::A1, plan::PlanT ; shift_factor::Real=one(T1)
	) where { T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2}, PlanT<:AbstractCCFPlan }
	num_lines_in_mask = length(plan.line_list)
	@assert num_lines_in_mask >= 1
	@assert size(projection,1) == length(λs)
	@assert size(projection,2) == 1 # Keep until implement CCF w/ multiple masks simultaneously
	projection .= zero(eltype(projection)) # Important to zero out projection before accumulating there!
	# set indexing variables
	p = 1
	p_left_edge_of_current_line = 1
	λsle_cur = λs[p]-0.5*(λs[p+1]-λs[p]) # Current pixel's left edge
	λsre_cur = 0.5*(λs[p]+λs[p+1]) # Current pixel's left edge
	#@assert issorted( plan.line_list.λ ) # Moved outside of inner loop
	m = searchsortedfirst(plan.line_list.λ, λsle_cur/shift_factor)
	while m <= length(plan.line_list.λ) && λ_min(plan.mask_shape,plan.line_list.λ[m]) * shift_factor < λsle_cur # only includemask entry if it fit entirely in chunk
	m += 1
	end
	if m > length(plan.line_list.λ) # Return early if no lines fall within chunk's lambda's
	return projection
	else
	#println("# Starting with mask entry at λ=", plan.line_list.λ[m], " for chunk with λ= ",first(λs)," - ",last(λs))
	end
	on_mask = false
	mask_mid = plan.line_list.λ[m] * shift_factor
	mask_hi = λ_max(plan.mask_shape,plan.line_list.λ[m]) * shift_factor # current line's upper limit
	mask_lo = λ_min(plan.mask_shape,plan.line_list.λ[m]) * shift_factor # current line's lower limit
	c_mps = speed_of_light_mps
	mask_weight = plan.line_list.weight[m]
	# loop through lines in mask, weight mask by fraction of PSF falling in each wavelength bin and divide by Δlnλ for pixel
	while m <= num_lines_in_mask
	if !on_mask
	if λsre_cur > mask_lo
	if λsre_cur > mask_hi # Pixel overshot this mask entry, so weight based on mask filling only a portion of the pixel,
	#For Tophat: projection[p,1] += mask_weight * (mask_hi - mask_lo) / (λsre_cur - λsle_cur)
	#frac_of_psf_in_v_pixel = 1.0
	one_over_delta_z_pixel = 0.5*(λsre_cur + λsle_cur) / (λsre_cur - λsle_cur)
	projection[p,1] += one_over_delta_z_pixel * mask_weight
	m += 1
	if m<=length(plan.line_list) # Move to next line
	mask_mid = plan.line_list.λ[m] * shift_factor
	mask_lo = λ_min(plan.mask_shape,plan.line_list.λ[m]) * shift_factor
	mask_hi = λ_max(plan.mask_shape,plan.line_list.λ[m]) * shift_factor
	mask_weight = plan.line_list.weight[m]
	else # We're out of lines, so exit
	break
	end
	else # Right edge of current pixel has entered the mask for this line, but left edge hasn't
	#For Tophat: projection[p,1] += mask_weight * (λsre_cur - mask_lo) / (λsre_cur - λsle_cur)
	frac_of_psf_in_v_pixel = integrate(plan.mask_shape, (mask_lo-mask_mid)/mask_mid*c_mps, (λsre_cur-mask_mid)/mask_mid*c_mps)
	one_over_delta_z_pixel = 0.5*(λsre_cur + λsle_cur) / (λsre_cur - λsle_cur)
	projection[p,1] += frac_of_psf_in_v_pixel * one_over_delta_z_pixel * mask_weight
	on_mask = true # Indicate next pixel will hav something to contribute based on the current line
	p_left_edge_of_current_line = p # Mark the starting pixel of the line in case the lines overlap
	p+=1 # Move to next pixel
	λsle_cur = λsre_cur # Current pixel's left edge
	λsre_cur = p<size(projection,1) ? 0.5*(λs[p]+λs[p+1]) : λs[p]+0.5*(λs[p]-λs[p-1]) # Current pixel's right edge
	end
	else # ALl of pixel is still to left of beginning of mask for this line.
	p+=1 # TODO: Opt try to accelerate with something like searchsortedfirst? Probably not worth it.
	λsle_cur = λsre_cur # Current pixel's left edge
	λsre_cur = p<size(projection,1) ? 0.5*(λs[p]+λs[p+1]) : λs[p]+0.5*(λs[p]-λs[p-1]) # Current pixel's right edge
	end
	else # on_mask == true
	if λsre_cur > mask_hi # Right edge of this pixel moved past the edge of the current line's mask
	#For Tophat: projection[p,1] += mask_weight * (mask_hi - λsle_cur) / (λsre_cur - λsle_cur)
	frac_of_psf_in_v_pixel = integrate(plan.mask_shape, (λsle_cur-mask_mid)*c_mps/mask_mid, (mask_hi-mask_mid)*c_mps/mask_mid)
	one_over_delta_z_pixel = 0.5*(λsre_cur + λsle_cur) / (λsre_cur - λsle_cur)
	projection[p,1] += frac_of_psf_in_v_pixel * one_over_delta_z_pixel * mask_weight
	on_mask = false # Indicate that we're done with this line
	m += 1 # Move to next line
	if m<=length(plan.line_list)
	mask_mid = plan.line_list.λ[m] * shift_factor
	mask_lo = λ_min(plan.mask_shape,plan.line_list.λ[m]) * shift_factor
	mask_hi = λ_max(plan.mask_shape,plan.line_list.λ[m]) * shift_factor
	mask_weight = plan.line_list.weight[m]
	if mask_lo < λsle_cur # If the lines overlap, we may have moved past the left edge of the new line. In that case go back to the left edge of the previous line.
	p = p_left_edge_of_current_line
	λsle_cur = p>1 ? 0.5*(λs[p]+λs[p-1]) : λs[p]-0.5*(λs[p+1]-λs[p]) # Current pixel's left edge
	λsre_cur = 0.5*(λs[p]+λs[p+1]) # Current pixel's right edge
	end
	else # We're done with all lines, can return early
	break
	end
	else # This pixel is entirely within the mask
	# assumes plan.mask_shape is normalized to integrate to unity and flux is constant within pixel
	#For Tophat: projection[p,1] += mask_weight # Add the full mask weight
	frac_of_psf_in_v_pixel = integrate(plan.mask_shape, (λsle_cur-mask_mid)/mask_mid*c_mps, (λsre_cur-mask_mid)/mask_mid*c_mps)
	one_over_delta_z_pixel = 0.5*(λsre_cur + λsle_cur) / (λsre_cur - λsle_cur)
	projection[p,1] += frac_of_psf_in_v_pixel * one_over_delta_z_pixel * mask_weight
	p += 1 # move to next pixel
	λsle_cur = λsre_cur # Current pixel's left edge
	λsre_cur = p<size(projection,1) ? 0.5*(λs[p]+λs[p+1]) : λs[p]+0.5*(λs[p]-λs[p-1]) # Current pixel's right edge
	end
	end
	if p>=size(projection,1) break end # We're done with all pixels in this chunk.
	end
	projection ./= c_mps
	return projection
	end