fist commit

.ipynb_checkpoints/Project-checkpoint.toml

@@ -0,0 +1,22 @@
+name = "Falcons"
+uuid = "eba0fb4b-2daf-4740-90d2-e5655143e0a0"
+authors = ["Yusuke Takase"]
+version = "0.1.0"
+
+[compat]
+julia = "1.5"
+Healpix = "2.3.0"
+StaticArrays = "1.2.0"
+ReferenceFrameRotations = "0.5.7"
+
+[deps]
+Healpix = "9f4e344d-96bc-545a-84a3-ae6b9e1b672b"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+ReferenceFrameRotations = "74f56ac7-18b3-5285-802d-d4bd4f104033"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]

Project.toml

@@ -4,7 +4,16 @@ authors = ["Yusuke Takase"]
 version = "0.1.0"
 
 [compat]
-julia = "1"
+julia = "1.5"
+Healpix = "2.3.0"
+StaticArrays = "1.2.0"
+ReferenceFrameRotations = "0.5.7"
+
+[deps]
+Healpix = "9f4e344d-96bc-545a-84a3-ae6b9e1b672b"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+ReferenceFrameRotations = "74f56ac7-18b3-5285-802d-d4bd4f104033"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/.ipynb_checkpoints/Falcons-checkpoint.jl

@@ -0,0 +1,16 @@
+module Falcons
+
+using ReferenceFrameRotations
+using Healpix
+using LinearAlgebra
+using Base.Threads
+using StaticArrays
+
+include("./function/func4scan.jl")
+include("./function/scanning.jl")
+include("./function/mapmake.jl")
+
+export get_scan_tod, scan_strategy
+export Mapmaking, Hitmap
+
+end

src/Falcons.jl

@@ -1,5 +1,16 @@
 module Falcons
 
-# Write your package code here.
+using ReferenceFrameRotations
+using Healpix
+using LinearAlgebra
+using Base.Threads
+using StaticArrays
+
+include("./function/func4scan.jl")
+include("./function/scanning.jl")
+include("./function/mapmake.jl")
+
+export get_scan_tod, scan_strategy
+export Mapmaking, Hitmap
 
 end

src/function/.ipynb_checkpoints/func4scan-checkpoint.jl

@@ -0,0 +1,16 @@
+function vec2ang_ver2(x, y, z)
+    norm = sqrt(x^2 + y^2 + z^2)
+    theta = acos(z / norm)
+    phi = atan(y, x)
+    if phi > π
+        phi = π - phi
+    end
+    return (theta, phi)
+end
+
+function quaternion_rotator(omega, t, rotate_axis)
+    #= Generate a quaternion that rotates by the angle omega*t around the rotate_axis axis. =#
+    rot_ang = @views omega * t
+    V = @views @SVector [cos(rot_ang/2.0), rotate_axis[1]*sin(rot_ang/2.0), rotate_axis[2]*sin(rot_ang/2.0), rotate_axis[3]*sin(rot_ang/2.0)]
+    return Quaternion(V)
+end

src/function/.ipynb_checkpoints/mapmake-checkpoint.jl

@@ -0,0 +1,55 @@
+function HitMap(nside, pix)
+    npix = nside2npix(nside)
+    hit_map = zeros(Int64, npix)
+    for i = eachindex(pix)
+        hit_map[pix[i]] += 1
+    end
+    return hit_map
+end
+
+
+function Mapmaking(scan_strategy_struct)
+    npix = nside2npix(scan_strategy_struct.nside)
+    split_num = 2
+    month = Int(scan_strategy_struct.times / split_num)
+    hwp_revol_rate = 2.0 * π * (scan_strategy_struct.hwp_rpm/ 60.0)
+
+    hit_map = zeros(Int32, npix)
+    Cross = zeros(Float32, (2,4, npix))
+    BEGIN = 0
+
+    @inbounds @simd for i = 1:split_num
+        println("process=", i, "/", split_num)
+        END = i * month
+
+        pix_tod, psi_tod = get_scan_tod(scan_strategy_struct, BEGIN, END)
+        println("Start mapmaking!")
+
+        @inbounds @simd for k = eachindex(pix_tod)
+            bore = @views pix_tod[k]
+            psi = @views psi_tod[k]
+            TIME = @views BEGIN + (k - 1) / scan_strategy_struct.sampling_rate
+
+            hit_map[bore] += 1
+
+            Cross[1,1,bore] += sin(psi)
+            Cross[2,1,bore] += cos(psi)
+            Cross[1,2,bore] += sin(2psi)
+            Cross[2,2,bore] += cos(2psi)
+            Cross[1,3,bore] += sin(3psi)
+            Cross[2,3,bore] += cos(3psi)
+            Cross[1,4,bore] += sin(4psi)
+            Cross[2,4,bore] += cos(4psi)
+
+        end
+        BEGIN = END + 1
+    end
+
+    link1 = @views @. (Cross[1,1,:]/hit_map)^2 + (Cross[2,1,:]/hit_map)^2
+    link2 = @views @. (Cross[1,2,:]/hit_map)^2 + (Cross[2,2,:]/hit_map)^2
+    link3 = @views @. (Cross[1,3,:]/hit_map)^2 + (Cross[2,3,:]/hit_map)^2
+    link4 = @views @. (Cross[1,4,:]/hit_map)^2 + (Cross[2,4,:]/hit_map)^2
+    out_map = @views [hit_map, link1, link2, link3, link4]
+
+    return out_map
+end

src/function/.ipynb_checkpoints/scanning-checkpoint.jl

@@ -0,0 +1,109 @@
+mutable struct scan_strategy
+    nside::Int
+    times::Int
+    sampling_rate::Int
+    alpha
+    beta
+    prec_period
+    spin_rpm
+    hwp_rpm
+    FP_theta
+    FP_phi
+    start_point
+    scan_strategy() = new()
+end
+
+
+@inline function get_scan_tod(scan_strategy_struct, start, stop)
+    nside = scan_strategy_struct.nside
+    alpha = @views deg2rad(scan_strategy_struct.alpha)
+    beta = @views deg2rad(scan_strategy_struct.beta)
+    FP_theta = scan_strategy_struct.FP_theta
+    FP_phi = scan_strategy_struct.FP_phi
+    omega_spin = @views (2π / 60) * scan_strategy_struct.spin_rpm
+    omega_prec = @views (2π / 60) / scan_strategy_struct.prec_period
+
+    #= Compute the TOD of the hitpixel and the TOD of the detector orientation at a specified sampling rate from time start to stop. =#
+
+    smp_rate = scan_strategy_struct.sampling_rate
+    loop_times = @views ((stop - start) * smp_rate)  + smp_rate
+    pix_tod = @views zeros(Int32, loop_times, length(FP_theta))
+    psi_tod = @views zeros(Float32, loop_times, length(FP_theta))
+    resol = @views Resolution(nside)
+
+    antisun_axis = @views @SVector [1.0, 0.0, 0.0]
+    spin_axis = @views @SVector [cos(alpha), 0.0, sin(alpha)]
+    y_axis = @views @SVector [0.0, 1.0, 0.0]
+    z_axis = @views @SVector [0.0, 0.0, 1.0]
+    omega_revol = @views (2π) / (60.0 * 60.0 * 24.0 * 365.25)
+
+    if scan_strategy_struct.start_point =="pole"
+        #= Set the initial position of the boresight to near the zenith. =#
+        bore_0 = @views @SVector [sin(alpha-beta), 0, cos(alpha-beta)]
+    end
+    if scan_strategy_struct.start_point == "equator"
+        #= Set the initial position of the boresight to near the equator. =#
+        bore_0 = @views @SVector [cos(alpha-beta), 0 , sin(alpha-beta)]
+    end
+    #= Generate the quaternion at the initial position of the boresight and detector orientation. =#
+    boresight_0 = @views Quaternion(0.0, bore_0)
+    detector_vec_0 = @views @SVector [-boresight_0.q3, 0.0, boresight_0.q1]
+    detector_orientation_0 = @views Quaternion(0.0, detector_vec_0)
+    travel_direction_vec_0 = @views bore_0 × detector_vec_0
+    travel_direction_0 = @views Quaternion(0.0, travel_direction_vec_0)
+
+    @views @inbounds @simd for j in eachindex(FP_theta)
+        #= 
+        Generating the pointhing quaternion of other detectors in the focal plane 
+        by adding an angular offset to the boresight based on the FP_theta and FP_phi information. 
+        =#
+        q_theta_in_FP = quaternion_rotator(deg2rad(FP_theta[j]), 1, y_axis)
+        q_phi_in_FP = quaternion_rotator(deg2rad(FP_phi[j]), 1, bore_0)
+        q_for_FP = q_phi_in_FP * q_theta_in_FP
+        pointing = q_for_FP * boresight_0 / q_for_FP
+        @views @inbounds @threads for i = 1:loop_times
+            t = start + (i - 1) / smp_rate
+            #= Generate the quaternion of revolution, precession, and spin at time t. =#
+            q_revol = quaternion_rotator(omega_revol, t, z_axis)
+            q_prec = quaternion_rotator(omega_prec, t, antisun_axis)
+            q_spin = quaternion_rotator(omega_spin, t, spin_axis)
+            #= 
+            P has the value of the pointing vector (x,y,z) at time t as its imaginary part.
+            The direction that the detector orientation is facing is then calculated as det_vec_t.
+            =#
+            Q = q_revol * q_prec * q_spin
+            P = Q * pointing / Q
+            q_det = Q * detector_orientation_0 / Q
+            travel_direction = Q * travel_direction_0 / Q
+            pointing_t = @SVector [P.q1, P.q2, P.q3]
+            travel_direction_vec = [travel_direction.q1, travel_direction.q2, travel_direction.q3]
+            #=
+            The direction of movement can be calculated by the outer product of the pointhing vector and the detector orientaion vector.
+            The vector of meridians can be calculated by combining the outer product of pointing and z axis.
+            =#
+            longitude = pointing_t × (pointing_t × z_axis)
+            bore_ang = vec2ang_ver2(pointing_t[1], pointing_t[2], pointing_t[3])
+
+            #ang_tod[1, i, j] = bore_ang[1]
+            #ang_tod[2, i, j] = bore_ang[2]
+
+            pix_tod[i, j] = ang2pixRing(resol, bore_ang[1], bore_ang[2])
+            #pix_tod[i, j] = ang2pix(_m_, bore_ang[1], bore_ang[2])
+
+            #=
+            The vector standing perpendicular to the sphere is defined as the outer product of the meridian and the moving_direction. 
+            The direction of divergent_vec depends on whether the trajectory enters the sphere from the left or the right side of the meridian.
+            =#
+            divergent_vec = longitude × travel_direction_vec
+            cosK = dot(travel_direction_vec, longitude) / (norm(travel_direction_vec) * norm(longitude))
+            #=
+            An if statement to prevent errors due to rounding errors that may result in |cosK|>1
+            =#
+            if abs(cosK) > 1.0
+                cosK = sign(cosK)
+            end
+            psi_tod[i, j] = acos(cosK) * sign(divergent_vec[3]) * sign(pointing_t[3])
+        end
+    end
+    return pix_tod, psi_tod
+end

src/function/func4scan.jl

@@ -0,0 +1,16 @@
+function vec2ang_ver2(x, y, z)
+    norm = sqrt(x^2 + y^2 + z^2)
+    theta = acos(z / norm)
+    phi = atan(y, x)
+    if phi > π
+        phi = π - phi
+    end
+    return (theta, phi)
+end
+
+function quaternion_rotator(omega, t, rotate_axis)
+    #= Generate a quaternion that rotates by the angle omega*t around the rotate_axis axis. =#
+    rot_ang = @views omega * t
+    V = @views @SVector [cos(rot_ang/2.0), rotate_axis[1]*sin(rot_ang/2.0), rotate_axis[2]*sin(rot_ang/2.0), rotate_axis[3]*sin(rot_ang/2.0)]
+    return Quaternion(V)
+end

src/function/mapmake.jl

@@ -0,0 +1,55 @@
+function HitMap(nside, pix)
+    npix = nside2npix(nside)
+    hit_map = zeros(Int64, npix)
+    for i = eachindex(pix)
+        hit_map[pix[i]] += 1
+    end
+    return hit_map
+end
+
+
+function Mapmaking(scan_strategy_struct)
+    npix = nside2npix(scan_strategy_struct.nside)
+    split_num = 2
+    month = Int(scan_strategy_struct.times / split_num)
+    hwp_revol_rate = 2.0 * π * (scan_strategy_struct.hwp_rpm/ 60.0)
+
+    hit_map = zeros(Int32, npix)
+    Cross = zeros(Float32, (2,4, npix))
+    BEGIN = 0
+
+    @inbounds @simd for i = 1:split_num
+        println("process=", i, "/", split_num)
+        END = i * month
+
+        pix_tod, psi_tod = get_scan_tod(scan_strategy_struct, BEGIN, END)
+        println("Start mapmaking!")
+
+        @inbounds @simd for k = eachindex(pix_tod)
+            bore = @views pix_tod[k]
+            psi = @views psi_tod[k]
+            TIME = @views BEGIN + (k - 1) / scan_strategy_struct.sampling_rate
+
+            hit_map[bore] += 1
+
+            Cross[1,1,bore] += sin(psi)
+            Cross[2,1,bore] += cos(psi)
+            Cross[1,2,bore] += sin(2psi)
+            Cross[2,2,bore] += cos(2psi)
+            Cross[1,3,bore] += sin(3psi)
+            Cross[2,3,bore] += cos(3psi)
+            Cross[1,4,bore] += sin(4psi)
+            Cross[2,4,bore] += cos(4psi)
+
+        end
+        BEGIN = END + 1
+    end
+
+    link1 = @views @. (Cross[1,1,:]/hit_map)^2 + (Cross[2,1,:]/hit_map)^2
+    link2 = @views @. (Cross[1,2,:]/hit_map)^2 + (Cross[2,2,:]/hit_map)^2
+    link3 = @views @. (Cross[1,3,:]/hit_map)^2 + (Cross[2,3,:]/hit_map)^2
+    link4 = @views @. (Cross[1,4,:]/hit_map)^2 + (Cross[2,4,:]/hit_map)^2
+    out_map = @views [hit_map, link1, link2, link3, link4]
+
+    return out_map
+end