backscatter_angle_sample3

A sample for the slice of the angular response.

''' A sample for the slice of the angular response.
'''

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

from irfpy.cena.backscatter_angle import fs



def integrate(sza, theta1):
    ''' Integrate the fs over full phi and theta=(0, theta1)
    '''

    n = 1000
    dphi = 360. / n
    dtheta = theta1 / n
    philist = np.linspace(0, 360, n, endpoint=False) + dphi / 2.
    thetalist = np.linspace(0, theta1, n, endpoint=False) + dtheta / 2.

    total = 0.

    for theta in thetalist:

        f = fs(sza, philist, theta)
        fsinthetadthetadphi = f.sum() * (dtheta * np.pi / 180.) * np.sin(theta * np.pi / 180.) * (dphi * np.pi / 180.0)
        total += fsinthetadthetadphi

    return total
        
def average(sza, theta1):
    integ = integrate(sza, theta1)
    omega = 2 * np.pi * (1 - np.cos(theta1 * np.pi / 180.))

    return integ / omega


if __name__ == '__main__':

    # Orbital slice: sza_at_equator = 45
    szaeq = 45
    szaeqrad = szaeq * np.pi / 180.
    lats = np.linspace(-100, 100, 200)

    phis = np.linspace(0, 360, 360)


    fig = plt.figure()
    ax = fig.add_subplot(111)

    for lat in lats:
        latrad = lat * np.pi / 180.
        locszarad = np.arccos(np.cos(latrad) * np.cos(szaeqrad))
        locsza = locszarad * 180. / np.pi

        flux = np.cos(locszarad)

        fs0 = fs(locsza, phis, 0) * flux

        ax.plot(np.ones_like(phis) * lat, fs0, '.')