# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# +
from __future__ import print_function
import os
import pickle

from shapely.geometry import LineString

from zipfile import ZipFile
import xml.sax, xml.sax.handler
from pyproj import Proj, transform

import numpy as np
import time
import datetime
from functools import reduce

from sklearn import svm, datasets
from sklearn.neural_network import MLPClassifier
from sklearn.preprocessing import StandardScaler

import scipy.io
import scipy.integrate
from astropy.convolution import convolve, Box1DKernel

from matplotlib import pyplot as plt
import matplotlib.cm as cm
from matplotlib import rc
from matplotlib.colors import LogNorm
import matplotlib.patches as patches
from matplotlib.collections import PatchCollection, LineCollection
from matplotlib.colors import ListedColormap, BoundaryNorm
from matplotlib.path import Path as mpath
from scipy.optimize import least_squares

from scipy import misc
from scipy.optimize import fsolve
import glob
from skimage import filters

import pandas as pd
from pathlib import Path

import warnings
warnings.filterwarnings('ignore')

# %matplotlib inline
pd.options.display.max_columns = 999
pd.options.display.max_rows = 90

np.set_printoptions(threshold=np.nan)

plt.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
## for Palatino and other serif fonts use:
#rc('font',**{'family':'serif','serif':['Palatino']})
plt.rc('text', usetex=True)

###############
## COMPUTER ##
###############

laptop = False

# +
k   = 0.4
f   = 7.2921*10**-5*2*np.sin(32*np.pi/180)
z0  = 0.01
A0  = 2
B0  = 4.5
G   = 10

z   = np.logspace(-2,4,500)

Ro  = G/f/z0
Ro

# +
k   = 0.4
f   = 7.2921*10**-5*2*np.sin(32*np.pi/180)
z0  = 0.01
A0  = 2
B0  = 4.5
G   = 10

def drag(us0):
    global A0, B0, f, z0, k, G
    return k**2/((np.log(us0/f/z0)-A0)**2+B0**2)-us0**2/G**2

us0 = fsolve(drag, G)
print(us0)

#find ug0 somehow then plot a0 (the deflection as a contormap with G and z0 as the axes, 
#will have to find us0 each time by solving the implicit function in a nested loop of G and z0
#this will give you the equilibrium deflection for a steady ABL
# i suppose i really want us and z0 to be the independent variables though...
a0  = np.arctan(f*B0/(np.log(us0/f/z0)-A0))
a0

# +
z0u = 10**-2
z0d = 10**0
Gud   = 10

G    = Gud

z0   = z0u
us0u = fsolve(drag, Gud)
a0u  = np.arcsin(us0u*B0/k/Gud)

z0   = z0d
us0d = fsolve(drag, Gud)
a0d  = np.arcsin(us0d*B0/k/Gud)
Hs   = us0d/f/50

X    = np.linspace(0,6000,100)
a0ud = a0u + iblc0*X**iblex/Hs*(a0d-a0u)

# +
k   = 0.4
f   = 7.2921*10**-5*2*np.sin(32*np.pi/180)
A0  = 2
B0  = 4.5

Z0  = np.logspace(-4,0,90)
GW  = np.linspace(0,50,100)

Z0g, GWg = np.meshgrid(Z0,GW)

US0 = np.empty_like(Z0g)
AL0 = np.empty_like(Z0g)
HS0 = np.empty_like(Z0g)

for i in np.arange(0,np.shape(Z0)[0]):
    for j in np.arange(0,np.shape(GW)[0]):
        
        z0 = Z0[i]
        G  = GW[j]
        
        us0 = fsolve(drag, G)
        
        US0[j,i] = us0   
        AL0[j,i] = np.arcsin(us0*B0/k/G)
        HS0[j,i] = us0/f/50
        
fetch = 6000
iblex = 1.4
iblc0 = 0.1/100
IBL = np.empty_like(Z0g)
IBL[:] = iblc0*fetch**iblex
# -

# %store -r Xmid_um
# %store -r DD

# +
plt.rcParams['text.usetex'] = True #Let TeX do the typsetting
plt.rcParams['text.latex.preamble'] = [r'\usepackage{sansmath}', r'\sansmath'] #Force sans-serif math mode (for axes labels)
plt.rcParams['font.family'] = 'sans-serif' # ... for regular text
plt.rcParams['font.sans-serif'] = 'Helvetica' # Choose a nice font here

S = 50
LW = 0.8
fs = 12
SC = 100

v1 = 0
v2 = 20
v3 = 270
v4 = 300

lpu1 = 0
hpu1 = 25

fig = plt.gcf()

ax3 = plt.subplot(221)

sc1 = ax3.pcolormesh(Z0g, GWg, US0, cmap='binary',zorder=-1,vmin=0,vmax=2)
sc2 = ax3.contour(Z0g, GWg, 180/np.pi*AL0, cmap='magma_r',zorder=-1,vmin=15,vmax=45)

ax3.set_xscale('log')

ax3.set_xlabel('$z_{0}$ (m)', fontsize=fs)
ax3.set_ylabel('$U_{g}$ (m/s)', fontsize=fs)

plt.xticks(fontsize=fs)
plt.yticks(fontsize=fs)

cbar3 = fig.colorbar(sc1, ticks=[0,1,2],ax=ax3, orientation='horizontal')
cbar3.ax.xaxis.set_ticks_position('bottom')
cbar3.set_label('$u_{*}$ (m/s)',fontsize=fs)
cbar3.ax.set_xticklabels(['0','1','2'])
cbar3.ax.xaxis.set_label_position('bottom')
cbar3.ax.tick_params(labelsize=fs)


ax3 = plt.subplot(222)

sc1 = ax3.pcolormesh(Z0g, GWg, 180/np.pi*AL0, cmap='magma_r',zorder=-1,vmin=15,vmax=45)
sc2 = ax3.contour(Z0g, GWg, US0, cmap='binary',zorder=-1,vmin=0,vmax=2)

ax3.set_xscale('log')

ax3.set_xlabel('$z_{0}$ (m)', fontsize=fs)
ax3.set_ylabel('$U_{g}$ (m/s)', fontsize=fs)

plt.xticks(fontsize=fs)
plt.yticks(fontsize=fs)

cbar3 = fig.colorbar(sc1, ticks=[15,30,45], ax=ax3, orientation='horizontal')
cbar3.ax.xaxis.set_ticks_position('bottom')
cbar3.set_label('$\\alpha$ ($^{\\circ}$)',fontsize=fs)
cbar3.ax.set_xticklabels(['15','30','45'])
cbar3.ax.xaxis.set_label_position('bottom')
cbar3.ax.tick_params(labelsize=fs)

fig.subplots_adjust(hspace=0.2)
fig.subplots_adjust(wspace=0.35)

ax3 = plt.subplot(223)

sc1 = ax3.pcolormesh(Z0g, GWg, HS0/IBL, cmap='seismic',zorder=-1,vmin=0,vmax=2)
sc3 = ax3.contour(Z0g, GWg, 180/np.pi*AL0, cmap='magma_r',zorder=-1,vmin=15,vmax=45)

ax3.set_xscale('log')

ax3.set_xlabel('$z_{0}$ (m)', fontsize=fs)
ax3.set_ylabel('$U_{g}$ (m/s)', fontsize=fs)

plt.xticks(fontsize=fs)
plt.yticks(fontsize=fs)

cbar3 = fig.colorbar(sc1, ticks=[0,1,2], ax=ax3, orientation='horizontal')
cbar3.ax.xaxis.set_ticks_position('bottom')
cbar3.set_label('$H_{Ek}/H_{SL}$',fontsize=fs)
cbar3.ax.set_xticklabels(['0','1','2'])
cbar3.ax.xaxis.set_label_position('bottom')
cbar3.ax.tick_params(labelsize=fs)

ax3 = plt.subplot(224)

sc1 = ax3.plot(X, 180/np.pi*a0ud,c='k')
sc1 = ax3.scatter(Xmid_um[Xmid_um>0],DD[Xmid_um>0]+90-10+180/np.pi*a0ud[0],c='y',s=1)

ax3.set_ylabel('$\\alpha$ ($^{\\circ}$)', fontsize=fs)
ax3.set_xlabel('$x$ (km)', fontsize=fs)

plt.xlim(0,6000)
plt.ylim(-45,135)

plt.xticks([0,2000,4000,6000], ('0', '2', '4', '6'))
plt.yticks([-45,0,45,90,135], ('-45', '0', '45', '90', '135'))

plt.xticks(fontsize=fs)
plt.yticks(fontsize=fs)

cbar3 = fig.colorbar(sc3, ax=ax3, orientation='horizontal')
cbar3.ax.xaxis.set_ticks_position('bottom')
cbar3.ax.xaxis.set_label_position('bottom')
cbar3.ax.tick_params(labelsize=fs)`

fig.subplots_adjust(hspace=0.2)
fig.subplots_adjust(wspace=0.35)

fig.set_size_inches(8, 10, forward=True)
# %cd /home/andrew/Documents/
plt.savefig('some_name86.png', bbox_inches='tight',dpi=300)
plt.savefig('some_name86.pdf', bbox_inches='tight')
# -

180/np.pi*a0d

StarcoderdataJupyterScriptsDedupFiltered