import Get_MW
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
import numpy as np

N=10 #number of frequencies

#These values are starting positions for coronal CME radio observations
ParmIn=29*[0]    # input array
ParmIn[0] =8e19  # Area, cm^2
ParmIn[1] =5e9  # Depth, cm
ParmIn[2] =3e6   # T_0, K
ParmIn[3] =0.05  # \eps (not used in this example)
ParmIn[4] =6.0   # \kappa (not used in this example)
ParmIn[5] =16    # number of integration nodes
ParmIn[6] =0.1   # E_min, MeV
ParmIn[7] =10.0  # E_max, MeV
ParmIn[8] =1.0   # E_break, MeV (not used in this example)
ParmIn[9] =4.0   # \delta_1
ParmIn[10]=6.0   # \delta_2 (not used in this example)
ParmIn[11]=1e8   # n_0 - thermal electron density, cm^{-3}
ParmIn[12]=1e6   # n_b - nonthermal electron density, cm^{-3}
ParmIn[13]=5.0 # B - magnetic field, G
ParmIn[14]=60.0  # theta - the viewing angle, degrees
ParmIn[15]=8.e7   # starting frequency to calculate spectrum, Hz
ParmIn[16]=0.005  # logarithmic step in frequency
ParmIn[17]=6     # Index of distribution over energy (KAP is chosen)
ParmIn[18]=N     # Number of frequencies (specified above)
ParmIn[19]=3     # Index of distribution over pitch-angle (GLC is chosen)
ParmIn[20]=90.0  # loss-cone boundary, degrees
ParmIn[21]=0.0   # beam direction (degrees) in GAU and SGA (not used in this example)
ParmIn[22]=0.2   # \Delta\mu
ParmIn[23]=0.0   # a_4 in SGA (not used in this example)
ParmIn[25]=12.0  # f^C_cr
ParmIn[26]=12.0  # f^WH_cr
ParmIn[27]=1     # matching on
ParmIn[28]=1     # Q-optimization on

def init_frequency_grid(startfreq,endfreq,numfreq=N):
    Params = ParmIn
    Params[16]=np.log10(endfreq/startfreq)/numfreq
    Params[15]=startfreq*1.e6
    Params[18]=numfreq
    s=Get_MW.GET_MW(Params)   # calling the main function
    f=s[0]           # emission frequency (GHz)
    fmhz=[i*1000. for i in f]
    return fmhz
    

def gs_kappa_func(freqgrid, temp=ParmIn[2],dens=ParmIn[11],kappa=ParmIn[4],emax=ParmIn[7],numfreq=N):
    Params = ParmIn
    Params[2]=temp
    Params[4]=kappa
    Params[7]=emax
    Params[11]=dens
    Params[15]=freqgrid[0]/1.e6
    Params[17]=6
    if not numfreq:
        numfreq=len(freqgrid)
    Params[16]=np.log10(freqgrid[-1]/freqgrid[0])/numfreq
    ParmIn[18]=numfreq
    s=Get_MW.GET_MW(ParmIn)   # calling the main function


    I_O=s[1]         # observed (at the Earth) intensity, O-mode (sfu)
    k_O=s[2]         # exp(-tau), O-mode
    #I_X=s[3]         # observed (at the Earth) intensity, X-mode (sfu)
    #k_X=s[4]         # exp(-tau), X-mode
    return I_O
