# 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()