#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on 01.01.2016

@author: richter bettina

Edited on 01.11.2022
@author: stephanie mayer
@author: richter bettina

"""


import numpy as np
import datetime
import pandas as pd

def process_profile(prof):

    ts = sorted(prof['data'].keys())
    for dt in ts:
        try:
            pro=prof['data'][dt]
            rc=rc_flat(pro)
            prof['data'][dt]['critical_cut_length'] = rc
            
            pendepth = comp_pendepth(pro, prof['info']['slopeAngle'])
            prof['data'][dt]['penetration_depth'] = pendepth
        except:
            continue

    return prof
#%%
    ##########################################################################################################################
    ########### Add improved and validated critical cut length for flat field according to Richter et al., 2019 (https://doi.org/10.5194/tc-13-3353-2019) ##############
    ########### Original formulation of Gaume et al., 2017 (https://doi.org/10.5194/tc-11-217-2017) resulted in negative values for slope simulations #################
    ########### In Mayer et al., 2022 (https://doi.org/10.5194/tc-16-4593-2022) flat rc values from slope simulation were ranked highly for computing probality of being unstable #########
    ########### Also add penetrations depth ##################################################################################
    ##########################################################################################################################

def rc_flat(pro):
    rho_ice = 917. #kg m-3
    gs_0 = 0.00125 #m
    thick = np.diff(np.concatenate((np.array([0]), pro['height'])))
    rho = pro['density']
    gs = pro['grain_size']*0.001 #[m]
    rho_sl = np.append(np.flip(np.cumsum(rho[::-1]*thick[::-1])/np.cumsum(thick[::-1]))[1:len(rho)],np.nan)
    #rho_sl_old = np.array([ np.sum( rho[i+1:]* thick[i+1:])/np.sum(thick[i+1:]) if i<len(rho)-1 else np.nan for i in range(len(rho)) ] )
    tau_p = pro['shear_strength']*1000. #[Pa]
    eprime = 5.07e9*(rho_sl/rho_ice)**5.13 / (1-0.2**2) #Eprime = E' = E/(1-nu**2) ; poisson ratio nu=0.2
    dsl_over_sigman = 1. / (9.81 * rho_sl) #D_sl/sigma_n = D_sl / (rho_sl*9.81*D_sl) = 1/(9.81*rho_sl)
    a = 4.6e-9
    b = -2.
    rc_flat = np.sqrt(a*( rho/rho_ice * gs/gs_0 )**b)*np.sqrt(2*tau_p*eprime*dsl_over_sigman)
    return rc_flat
        
def comp_pendepth(pro, slopeangle):
    top_crust = 0
    thick_crust = 0
    rho_Pk = 0
    dz_Pk = 1.e-12
    crust = False
    e_crust = -999
    
    layer_top = pro['height']
    ee = len(layer_top)-1
    thick = np.diff(np.concatenate((np.array([0]), layer_top)))
    rho = pro['density']
    HS = layer_top[-1]
    graintype = pro['grain_type']  
    min_thick_crust = 3 #cm
    
    while (ee >= 0) & ((HS-layer_top[ee])<30):
        
        rho_Pk = rho_Pk + rho[ee]*thick[ee]
        dz_Pk = dz_Pk + thick[ee]
        
        if crust == False:
        ##Test for strong mf-crusts MFcr.
        ## Look for the first (from top) with thickness perp to slope > 3cm
            if (graintype[ee] == 772) & (rho[ee] >500.): ## very high density threshold, but implemented as this in SP
                if e_crust == -999:
                   e_crust = ee
                   top_crust = layer_top[ee]
                   thick_crust = thick_crust + thick[ee]
                elif (e_crust - ee) <2:
                   thick_crust = thick_crust + thick[ee]
                   e_crust = ee
            elif e_crust > 0:
               if thick_crust*np.cos(np.deg2rad(slopeangle)) > min_thick_crust:
                   crust = True
               else:
                   e_crust = -999
                   top_crust = 0
                   thick_crust = 0

        ee = ee-1
                         

    
    rho_Pk = rho_Pk/dz_Pk        #average density of the upper 30 cm slab
    return np.min([0.8*43.3/rho_Pk, (HS-top_crust)/100.]) #NOTE  Pre-factor 0.8 introduced May 2006 by S. Bellaire , Pk = 34.6/rho_30  # factor "/100" to get upper crust depth in meter
