gridsphere_spacecraft_fov_surface

Spacecraft-Surface field of view mapping sample

''' Spacecraft-Surface field of view mapping sample
'''
import os
import sys
import logging
logging.basicConfig()
import datetime
import math

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

import irfpy.util.gridsphere


def spacecraft_position():
    n = 100
    rd = 1.1  # 10% bigger than the sphere

    angle = np.linspace(0, 360, n) * np.pi / 180.0

    x = np.cos(angle) * np.cos(25.8 * np.pi / 180.0) * rd
    y = np.cos(angle) * np.sin(25.8 * np.pi / 180.0) * rd
    z = np.sin(angle) * rd

    tilt = np.array([[1, 0, 0],
                    [0, np.cos(0.03), np.sin(0.03)],
                    [0, -np.sin(0.03), np.cos(0.03)]])

    xyz = np.array([x, y, z])    # 3, n
    xyz = tilt.dot(xyz)

    x, y, z = xyz

    r = np.empty_like(x)
    r.fill(rd)
    lat = np.arcsin(z / rd) * 180. / np.pi
    lon = np.arctan2(y, x) * 180. / np.pi

    while lon.max() >= 360:
        lon = np.where(lon >= 360, lon - 360, lon)
    while lon.min() < 0:
        lon = np.where(lon < 0, lon + 360, lon)

    return (np.array([x, y, z]), np.array([r, lon, lat]))


def main1():

    fov = 20   # 20 degrees.
    cosfov = np.cos(fov * np.pi / 180.)

    xyz, geo = spacecraft_position()

    # Simple implementation, 180x360
    sp = irfpy.util.gridsphere.SimpleGridSphere(nlon=360, nlat=180)

    
    # xyz, r, lon, lat is the spacecraft position.
    for (x, y, z), (r, lon, lat) in zip(xyz.T, geo.T):
        print((x, y, z), (r, lon, lat))

        xyz = np.array([x, y, z])  # Position vector
        xyz_norm = xyz / r

        visible_angle = np.arcsin(r * np.sin(fov * np.pi / 180.)) - fov * np.pi / 180.0
        cos_visible = np.cos(visible_angle)

        ngrid = len(sp)   # Sphere element size.

        for igrid in range(ngrid):
            # Loop for each position on the sphere.
            slon, slat = sp.index2lonlat(igrid)

            slon = slon * (np.pi / 180.0)
            slat = slat * (np.pi / 180.0)

            # Surface element position in a vector form
            xyz_surf = sp.index2xyz(igrid)

            # Visibility.
            if xyz_surf.dot(xyz_norm) < cos_visible:
                continue

            # View vector of the surface element from spacecraft
            viewvec_surf_sc = xyz_surf - xyz
            viewvec_surf_sc = viewvec_surf_sc / np.linalg.norm(viewvec_surf_sc)

            cosine = viewvec_surf_sc.dot(-xyz_norm)

            if cosine < cosfov:
                continue

            print(igrid, slon * 180 / np.pi, slat * 180 / np.pi)
            sp.append_value(igrid, 1)

    
    lons, lats = sp.get_bgrid()
    ave = sp.get_average()
    cnt = sp.get_counts()

    plt.pcolor(lons, lats, cnt)
    plt.plot(geo[1, :], geo[2, :], 'r+')


    plt.savefig('test.png')


if __name__ == '__main__':
    main1()