#!/usr/bin/env python3
# coding: utf-8
"""
  NAME: 
    AerosolSpectrum-DavidDelene.py  

  PURPOSE:
    Open .csv files containing CDP and PCASP Channel Number, Minimum Diameter, 
    Maximum Diameter, and Concentration data to calculate the number concentration, 
    area distribution, and volume distribution curves. The script also calculates the 
    total water content, total PM2.5 mass, number concentration, and area to be included 
    in matplotlib bar plots for each distribution

  SYNTAX: 
    Input: Program REQUIRES input CDP and PCASP .csv files to run. 
    
    Output: The program will output 6 total .png files. Three .png files 
    will be the number concentration, area distribution, volume 
    distribution for the CDP while the following three .png images will be for the PCASP data. 

  PARAMETERS:
    *None* 

  EXECUTION EXAMPLE:
    python ATSC565_Air_Quality_HW1_AS.py

  MODIFICATION HISTORY:
    David Delene - David Delene <david.delene.und@gmail.com> - 2025/11/05
      Written.

  Copyright 2025 David Delene

  This program is distributed under the terms of the GNU General Public License

  This file is part of Airborne Data Processing and Analysis (ADPAA).

  ADPAA is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  ADPAA is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with ADPAA.  If not, see <http://www.gnu.org/licenses/>.
"""

import numpy as np
import pandas as pd 
import sys
import matplotlib.pyplot as plt

# === Calculate Aerosols Spectrum. ===
cdp_data = pd.read_csv('./CDP_data.csv', sep=',')

rho_water = 1e6  # g/m3

channel_number_cdp = cdp_data['Channel Number'].tolist()        # Channel number
min_diameter_cdp   = cdp_data['Min Diameter (um)'].tolist()     # micron
max_diameter_cdp   = cdp_data['Max Diameter (um)'].tolist()     # micron
concentration_cdp  = cdp_data['Concentration (#/cm3)'].tolist() # #/cm3

mean_diameter_cdp        = []
number_concentration_cdp = []
area_distribution_cdp    = []
volume_distribution_cdp  = []
str_list                 = []
x_label                  = []
concentration_val        = []
concentration_val_m      = []
mean_diam                = []
mean_diam_m              = []
denominator_cdp          = []
number_concentration_cdp = []
area_diameter            = []
area_distribution        = []
volume_diameter          = []
volume_distribution      = []
lwc                      = []

for i in range(len(channel_number_cdp)):
    min_diam = float(min_diameter_cdp[i]) # micron
    max_diam = float(max_diameter_cdp[i]) # micron
    
    str_diam = str(min_diam)+'-'+str(max_diam)
    str_list.append(str_diam)
    x_label.append(str((max_diam+min_diam)/2.0))
    
    concentration_val.append(float(concentration_cdp[i])) # #/cm3
    concentration_val_m.append(float(concentration_cdp[i]*(1e6))) # #/m3

    mean_diam.append((float(max_diameter_cdp[i]) + float(min_diameter_cdp[i]))/2.0) # micron
    mean_diam_m.append(((float(max_diameter_cdp[i]) + float(min_diameter_cdp[i]))/2.0)*(1e-6)) # micron
    
    denominator_cdp.append(np.log(float(max_diameter_cdp[i])) - np.log(float(min_diameter_cdp[i])))
    
    number_concentration_cdp.append(concentration_val[i]/denominator_cdp[i]) # [dN/d log Dp [cm-3]]
    
    area_diameter.append(concentration_val[i] * np.pi * (mean_diam[i]**2))
    area_distribution_cdp.append(area_diameter[i]/denominator_cdp[i]) # dA/d log Dp [um2cm-3] 
    
    volume_diameter.append(concentration_val[i] * (np.pi/6) * (mean_diam[i]**3))  # [um3]
    volume_distribution_cdp.append(volume_diameter[i]/denominator_cdp[i])  # [um3 cm-3]

    lwc.append((np.pi/6) * 1e-6 * concentration_cdp[i] * (mean_diam[i])**3)

lwc_total_val = np.sum(lwc)
total_number_con = np.sum(concentration_cdp)
total_area_dist = np.sum(area_diameter)

fig=plt.figure()
plt.plot(x_label, number_concentration_cdp, zorder=2, color='blue')
plt.text(x_label[0], 900, 'Total Number Concentration: %s'%(str(total_number_con)[:6]))
plt.xlabel('Diameter (micron)')
plt.ylabel('Number Concentration (dN/d log Dp '+'$[cm^{-3}])$')
plt.grid(alpha=0.35, zorder=1)
plt.savefig('cdp_number_concentration.png', bbox_inches='tight')
plt.close()

fig = plt.figure()
plt.plot(x_label, area_distribution_cdp, zorder=2, color='red')
plt.text(x_label[0], 850000, 'Total Particle Area: %s um$^2$'%(str(total_area_dist)[:6]))
plt.yscale('log')
plt.xlabel('Diameter (micron)')
plt.ylabel('$Area Distribution (dA/d log Dp [um^2 cm^{-3}])$')
plt.grid(alpha=0.35, zorder=1)
plt.savefig('cdp_area_distribution.png', bbox_inches='tight')
plt.close()

fig = plt.figure()
plt.plot(x_label, volume_distribution_cdp, zorder=2, color='black')
plt.text(x_label[0], 1500000, 'Total Liquid Water Content (g/m$^3$): %s'%(str(lwc_total_val)[:4]))
plt.yscale('log') 
plt.xlabel('Diameter (micron)')
plt.ylabel('Volume Distribution (dV/d log Dp [um$^3 cm${-3}])')
plt.grid(alpha=0.35, zorder=1)
plt.savefig('cdp_volume_distribution.png', bbox_inches='tight')
plt.close()

'''# === Calculate Aerosols Spectrum. ===
pcasp_data = pd.read_csv('./PCASP_data.csv', sep=',')

rho_ammon = 1.769  # g/cm3

channel_number_pcasp = pcasp_data['Channel Number'].tolist()  # Channel number
min_diameter_pcasp = pcasp_data['Min Diameter (um)'].tolist()  # micron
max_diameter_pcasp = pcasp_data['Max Diameter (um)'].tolist()  # micron
concentration_pcasp = pcasp_data['Concentration (#/cm3)'].tolist() # #/cm3

mean_diameter_pcasp = []
number_concentration_pcasp = []
area_distribution_pcasp = []
volume_distribution_pcasp = []

str_list_pcasp = []

for i in range(len(channel_number_pcasp)):
    
    max_diam = float(max_diameter_pcasp[i]) # micron
    min_diam = float(min_diameter_pcasp[i]) # micron
    
    str_diam = str(min_diam)+'-'+str(max_diam)
    str_list_pcasp.append(str_diam)
    
    concentration_val = float(concentration_pcasp[i]) # #/cm3
    concentration_val_m = concentration_val*(1e6) # #/m3
    
    mean_diam = (max_diam + min_diam)/2 # micron
    mean_diam_m = mean_diam * (1e-6) # meters
    mean_diameter_pcasp.append(mean_diam)
    
    denominator_pcasp = np.log(max_diam) - np.log(min_diam)
    
    number_concentration_val = concentration_val/denominator_pcasp  #[dN/d log Dp [cm-3]]
    number_concentration_pcasp.append(number_concentration_val)
    
    area_diameter = concentration_val * np.pi * (mean_diam**2)  #dA
    area_dist_val = area_diameter/denominator_pcasp   #dA/d logDp [um2cm-3]
    area_distribution_pcasp.append(area_dist_val)
    
    volume_diameter = concentration_val * (np.pi/6) * (mean_diam**3)  # [um3]
    vol_dist_val = volume_diameter/denominator_pcasp  # [um3 cm-3]
    volume_distribution_pcasp.append(vol_dist_val)

pm_total_val = rho_ammon * (np.pi/6) *np.sum(np.array(concentration_pcasp) * (np.array(mean_diameter_pcasp)**3))
total_number_con = sum(number_concentration_pcasp)
total_area_dist = sum(area_distribution_pcasp)

fig=plt.figure()
plt.bar(str_list_pcasp, number_concentration_pcasp, zorder=2, color='blue')
plt.text(str_list_pcasp[-5], 900, 'Total Number \nConcentration: %s'%(str(total_number_con)[:6]))
plt.yscale('log')
plt.xlabel('Diameter (micron)')
plt.ylabel('Number Concentration (dN/d log Dp [cm-3])')
plt.xticks(rotation=90)
plt.grid(alpha=0.35, zorder=1)
plt.savefig('pcasp_number_concentration.png', bbox_inches='tight')
plt.close()

fig = plt.figure()
plt.bar(str_list_pcasp, area_distribution_pcasp, zorder=2, color='red')
plt.text(str_list_pcasp[-5], 75, 'Total Particle \nArea: %s'%(str(total_area_dist)[:9]))
plt.yscale('log')
plt.xlabel('Diameter (micron)')
plt.ylabel('Area Distribution (dA/d log Dp [um2 cm-3])')
plt.xticks(rotation=90)
plt.grid(alpha=0.35, zorder=1)
plt.savefig('pcasp_area_distribution.png', bbox_inches='tight')
plt.close()

fig = plt.figure()
plt.bar(str_list_pcasp, volume_distribution_pcasp, zorder=2, color='black')
plt.text(str_list_pcasp[0], 12, 'Total PM2.5 \nMass (ug/m3): %s'%(str(pm_total_val)[:5]))
plt.yscale('log')
plt.xlabel('Diameter (micron)')
plt.ylabel('Volume Distribution (dV/d log Dp [um3 cm-3])')
plt.xticks(rotation=90)
plt.grid(alpha=0.35, zorder=1)
plt.savefig('pcasp_volume_distribution.png', bbox_inches='tight')
plt.close()
'''