backscatter_angle_sample2

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 average2(sza, theta_range, phi_range, n=50):
    ''' Integrate and average the fs over theta_range and phi_range, in degrees.
    '''
    p0, p1 = phi_range
    t0, t1 = theta_range
    dphi = (p1 - p0) / n
    dtheta = (t1 - t0) / n

    philist = np.linspace(p0, p1, n, endpoint=False) + dphi / 2.
    thetalist = np.linspace(t0, t1, n, endpoint=False) + dtheta / 2.

    total_f_str = 0.
    total_str = 0.

    for theta in thetalist:
        f = fs(sza, philist, theta)
        total_f_str += f.sum() * (dtheta * np.pi / 180.0) * np.sin(theta * np.pi / 180.0) * (dphi * np.pi / 180.0)
        total_str += n * (dtheta * np.pi / 180.0) * np.sin(theta * np.pi / 180.0) * (dphi * np.pi / 180.0)

    return total_f_str / total_str


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__':

    # should be the same.
    av15 = average(45.0, 22.5)
    av15_2 = average2(45.0, [0, 22.5], [0, 360.])
    print(av15, av15_2)

    fig = plt.figure()
    ax = fig.add_subplot(111)
    for sza in np.arange(0., 91, 10):
        f90 = average2(sza, [0, 90.], [0, 360.])
        for phi in np.arange(0., 360, 5):
            f15 = average2(sza, [0, 22.5], [phi, phi+1])
            ax.plot([sza], [f90 / f15 * 2 * np.pi], 'b.')

    szas = np.arange(0, 91, 1)
    ratios = np.zeros_like(szas) * 1.0
    for i, sza in enumerate(szas):
        f15 = integrate(sza, 22.5)
        av15 = average(sza, 22.5)
        f90 = integrate(sza, 90.0)
        av90 = average(sza, 90.0)
#        ax.plot([sza], [av90/ av15 * 2 * np.pi], 'r+')
        ratios[i] = av90 / av15 * 2 * np.pi
        print(sza, f15, f90, f90 / f15, av15, av90, av90 / av15, av90 / av15 * 2 * np.pi, ratios[i])

    print(szas, ratios)
    ax.plot(szas, ratios, 'r')

    ax.plot([0, 90], [2 * np.pi, 2 * np.pi], 'g')

    ax.set_ylim(0, 15)

    fig.savefig('backscatter_angle_sample2.eps')