<jupyter_start><jupyter_text>Samara Plots NotebookI just wanted to pull together some quick plots of the data into one place, so I'll do so inline here.First, I'll write down where the data is in this repository for my benefit. Directory layout (data locations)* Density Study    * Samara Segment Mass.xlsx* DropStudy    * OutdoorDropData        * rawdataworkspace.mat        * 2017-04 Samara Drop Tests.xlsx        * weatherData.xlsx    * VideoScripts        * NorwayMapleData.csv        * RealMapleData.csv        * SilverMapleData.csv         Load data:<jupyter_code>real_vel = pd.read_csv('DropStudy/VideoScripts/RealMapleData.csv',header=0,names=['trial','cutoff','velocity','rotation'])
norway_vel = pd.read_csv('DropStudy/VideoScripts/NorwayMapleData.csv',header=0,names=['trial','cutoff','velocity','rotation'])
silver_vel = pd.read_csv('DropStudy/VideoScripts/SilverMapleData.csv',header=0,names=['trial','cutoff','velocity','rotation'])

vels = list(norway_vel['velocity'].values)
angs = list(norway_vel['rotation'].values)
Nvelocities = []; 
Nrotation  = []; 
for val in zip(vels,angs):
    if not ' ' in val:
        Nvelocities.append(float(val[0]))
        Nrotation.append(float(val[1]))

print('3D-printed Norway maples \nMean (Vel): ',np.mean(Nvelocities),'\n','Std (Vel): ',np.std(Nvelocities),'\n''Mean (Rot): ',np.mean(Nrotation),'\n','Std (Rot): ',np.std(Nrotation),'\n\n')
plt.figure()
plt.hist(Nvelocities)
plt.figure()
plt.hist(Nrotation)
plt.figure(3, figsize=(3.2,2.8))
plt.scatter(Nrotation,Nvelocities, marker='^', edgecolors='k', s=16, c=np.array([[0,1,0]]))
plt.errorbar(np.mean(Nrotation), np.mean(Nvelocities), xerr=np.std(Nrotation), yerr=np.std(Nvelocities), fmt='^', color='black', ms=6)

vels =list(silver_vel['velocity'].values)
angs =list(silver_vel['rotation'].values)
Svelocities  = []; 
Srotation = []; 
for val in zip(vels,angs):
    if not ' ' in val:
        Svelocities.append(float(val[0]))
        Srotation.append(float(val[1]))

print('3D-printed Silver maples \nMean (Vel): ',np.mean(Svelocities),'\n','Std (Vel): ',np.std(Svelocities),'\n''Mean (Rot): ',np.mean(Srotation),'\n','Std (Rot): ',np.std(Srotation),'\n\n')

plt.figure()
plt.hist(Svelocities)
plt.figure()
plt.hist(Srotation)
plt.figure(3)
plt.scatter(Srotation,Svelocities, marker = 'v', edgecolors='k', s=16, c=np.array([[0.8,0.8,0]]))
plt.errorbar(np.mean(Srotation), np.mean(Svelocities), xerr=np.std(Srotation), yerr=np.std(Svelocities), fmt='v', color='black', ms=6)

vels =list(real_vel['velocity'].values)
angs =list(real_vel['rotation'].values)
Rvelocities  = []; 
Rrotation  = []; 
for val in zip(vels,angs):
    if not ' ' in val:
        Rvelocities.append(float(val[0]))
        Rrotation.append(float(val[1]))

print('Real Norway maples \nMean (Vel): ',np.mean(Rvelocities),'\n','Std (Vel): ',np.std(Rvelocities),'\n''Mean (Rot): ',np.mean(Rrotation),'\n','Std (Rot): ',np.std(Rrotation),'\n\n')

plt.figure()
plt.hist(Rvelocities)
plt.figure()
plt.hist(Rrotation)
plt.figure(3)
plt.scatter(Rrotation,Rvelocities, edgecolors='k', s=16, c=0.9*np.array([[1,165/255,0]]))
plt.errorbar(np.mean(Rrotation), np.mean(Rvelocities), xerr=np.std(Rrotation), yerr=np.std(Rvelocities), fmt='o', color='black', ms=6)
plt.xlabel('Rotational velocity (rad/sec)', fontsize=14)
plt.ylabel('Descent speed (m/s)', fontsize=14)
plt.gca().tick_params(labelsize=11)
plt.savefig('Figures/IndoorDropData.png', dpi=450, bbox_inches='tight', transparent=True)
plt.savefig('Figures/IndoorDropData.eps', bbox_inches='tight', transparent=True)
plt.show()

#plt.scatter(Nrotation,Nvelocities, marker='^', edgecolors='k', s=16, c=np.array([[0,1,0]]))
#plt.errorbar(np.mean(Nrotation), np.mean(Nvelocities), xerr=np.std(Nrotation), yerr=np.std(Nvelocities), fmt='^', color='black', ms=6)

#plt.scatter(Srotation,Svelocities, marker = 'v', edgecolors='k', s=16, c=np.array([[0.8,0.8,0]]))
#plt.errorbar(np.mean(Srotation), np.mean(Svelocities), xerr=np.std(Srotation), yerr=np.std(Svelocities), fmt='v', color='black', ms=6)

plt.scatter(Rrotation,Rvelocities, edgecolors='k', s=16, c=0.9*np.array([[1,165/255,0]]))
plt.errorbar(np.mean(Rrotation), np.mean(Rvelocities), xerr=np.std(Rrotation), yerr=np.std(Rvelocities), fmt='o', color='black', ms=6)

plt.xlabel('Rotational velocity (rad/sec)', fontsize=14)
plt.ylabel('Descent speed (m/s)', fontsize=14)
plt.gca().tick_params(labelsize=11)

import pandas as pd
import numpy as np

raw_od_data = pd.read_excel('DropStudy/OutdoorDropData/2017-04 Samara Drop Tests.xlsx', header=1, usecols=[a for a in range(13)])

new_df = pd.DataFrame(columns=['SamaraType', 'Group', 'Trial', 'Date', 'Time', 'Dist', 'Ang'])

samara_list = ['SM', 'NM', 'N']
prop_list = ['Time', 'Dist', 'Ang']

checked_list = []

for index, row in raw_od_data.iterrows():
#   print(row)
    for samara_type in samara_list:
        group = row['Drop']
        for trial in range(1,4):
            if (samara_type, group, trial) not in checked_list:
                checked_list.append((samara_type, group, trial))
                break
        row_dict = {'SamaraType': samara_type, 'Group': row['Drop'], 'Trial': trial, 'Date': row['Date']}
        for prop in prop_list:
            value = row[prop + ' - ' + samara_type]
            if not value == '-' and ~np.isnan(value):
                row_dict[prop] = value
        if len(row_dict) > 3:
            new_df = new_df.append(pd.Series(row_dict), ignore_index=True)
        else:
            checked_list = checked_list[:-1]
#         new_df.append(pd.DataFrame(row_dict), ignore_index=True)
    #         print(index, prop, samara_type, val)
new_df['Velocity'] = 10*0.3048/new_df['Time']
fig1, ax1 = plt.subplots(figsize=(3.2,2.8))
fig2, ax2 = plt.subplots(nrows=1, ncols=3, figsize=(6.5,2.8), sharex=True, sharey=True)

for row in new_df.itertuples():
    if ~np.isnan(row.Velocity):
        if row.SamaraType == 'N':
            color = 0.9*np.array([[1,165/255,0]])
            marker = 'o'
        elif row.SamaraType == 'NM':
            color = np.array([[0,1,0]])
            marker = '^'
        else:
            color = np.array([[0.8,0.8,0]])
            marker = 'v'
        ax1.scatter(row.Velocity, row.Dist, edgecolors='k', s=16, marker=marker, c=color)
        date = row.Date
        if date==pd.Timestamp('2017-04-24'):
            ax2[0].scatter(row.Dist*np.sin(row.Ang*np.pi/180)*0.0254, row.Dist*np.cos(row.Ang*np.pi/180)*0.0254, marker=marker, edgecolors='k', s=16, c=color)
        elif date==pd.Timestamp('2017-04-26'):
            ax2[1].scatter(row.Dist*np.sin(row.Ang*np.pi/180)*0.0254, row.Dist*np.cos(row.Ang*np.pi/180)*0.0254, marker=marker, edgecolors='k', s=16, c=color)
        elif date==pd.Timestamp('2017-04-27'):
            ax2[2].scatter(row.Dist*np.sin(row.Ang*np.pi/180)*0.0254, row.Dist*np.cos(row.Ang*np.pi/180)*0.0254, marker=marker, edgecolors='k', s=16, c=color)
plt.sca(ax1)
plt.xlabel('Descent speed (m/s)', fontsize=14)
plt.ylabel('Distance traveled (m)', fontsize=14)
ax1.tick_params(labelsize=11)
# plt.savefig('Figures/OutdoorSpeedDist.png', dpi=450, bbox_inches='tight', transparent=True)
# plt.savefig('Figures/OutdoorSpeedDist.eps', bbox_inches='tight', transparent=True)
plt.sca(ax2[0])
plt.ylabel('Landing position (m)', fontsize=14)
plt.sca(ax2[1])
plt.xlabel('Landing position (m)', fontsize=14)
for k in range(3):
    ax2[k].tick_params(labelsize=11)
    ax2[k].set_aspect('equal')
# plt.savefig('Figures/OutdoorScatter.png', dpi=450, bbox_inches='tight', transparent=True)
# plt.savefig('Figures/OutdoorScatter.eps', bbox_inches='tight', transparent=True)
plt.show()
collected_data = pd.DataFrame(columns=['SamaraType','Group',
                                       'ODspeedMean','ODspeedStd','ODdistMean','ODdistStd',
                                       'IDspeedMean','IDspeedStd','IDrotMean','IDrotStd'])

fig1, ax1 = plt.subplots()
fig2, ax2 = plt.subplots()

for samType in ['SM', 'NM', 'N']:
    odvels = []
    if samType=='SM':
        marker = 'v'
        color = 'yellow'
        prefix = 's'
        IDdata = silver_vel.copy(deep=False)
    elif samType=='NM':
        marker = '^'
        color = 'green'
        prefix = 'n'
        IDdata = norway_vel.copy(deep=False)
    else:
        marker = 'o'
        color = 'orange'
        prefix = 'r'
        IDdata = real_vel.copy(deep=False)
    ODdata = new_df[new_df['SamaraType']==samType]
    for group in range(1,31):
        ODgroupData = ODdata[ODdata['Group']==group]
        temp_row = {'SamaraType':prefix,'Group':group}
        ODspeeds = []
        ODdist = []
        IDspeeds = []
        IDrotation = []
        for row in ODgroupData.itertuples():
            trial = row.Trial
            ODspeeds.append(row.Velocity)
            odvels.append(row.Velocity)
            ODdist.append(row.Dist)
        for trial in range(1,4):
            key = '{:s}-g{:02d}-t{:02d}-data.csv'.format(prefix,group,trial)
            if len(IDdata[IDdata['trial']==key])> 0:
                if not IDdata[IDdata['trial']==key]['velocity'].values[0]==' ':
                    IDspeeds.append(float(IDdata[IDdata['trial']==key]['velocity'].values[0]))
                if not IDdata[IDdata['trial']==key]['rotation'].values[0]==' ':
                    IDrotation.append(float(IDdata[IDdata['trial']==key]['rotation'].values[0]))
        if len(ODspeeds)>0:
            temp_row['ODspeedMean'] = np.nanmean(ODspeeds)
            temp_row['ODspeedStd'] = np.nanstd(ODspeeds)
            temp_row['ODdistMean'] = np.nanmean(ODdist)
            temp_row['ODdistStd'] = np.nanstd(ODdist)
        else:
            temp_row['ODspeedMean'] = np.nan
            temp_row['ODspeedStd'] = np.nan
            temp_row['ODdistMean'] = np.nan
            temp_row['ODdistStd'] = np.nan
        if len(IDspeeds)>0:
            temp_row['IDspeedMean'] = np.nanmean(IDspeeds)
            temp_row['IDspeedStd'] = np.nanstd(IDspeeds)
            temp_row['IDrotMean'] = np.nanmean(IDrotation)
            temp_row['IDrotStd'] = np.nanstd(IDrotation)
        else:
            temp_row['IDspeedMean'] = np.nan
            temp_row['IDspeedStd'] = np.nan
            temp_row['IDrotMean'] = np.nan
            temp_row['IDrotStd'] = np.nan
        collected_data = collected_data.append(pd.Series(temp_row), ignore_index=True)
        
        ax1.errorbar(temp_row['IDspeedMean'], temp_row['ODspeedMean'], 
                     xerr=temp_row['IDspeedStd'], yerr=temp_row['ODspeedStd'], 
                     marker=marker, color=color, ecolor='k', mec='k', ms=8)
        
        ax2.errorbar(temp_row['IDrotMean'], temp_row['ODdistMean']*0.0254, 
                     xerr=temp_row['IDrotStd'], yerr=temp_row['ODdistStd']*0.0254, 
                     marker=marker, color=color, ecolor='k', mec='k', ms=8)
    print(samType, np.nanmean(odvels), '$\pm$', np.nanstd(odvels))

ax1.set_aspect('equal')
plt.sca(ax1)
plt.xlim([0,2.8])
plt.ylim([0,2.8])
plt.xlabel('Indoor velocity (m/s)', fontsize=12)
plt.ylabel('Outdoor velocity (m/s)', fontsize=12)
plt.savefig('IOspeedcomp.eps', transparent=True, bbox_inches='tight')
plt.savefig('IOspeedcomp.png', dpi=300, transparent=True, bbox_inches='tight')


ax1.set_aspect('equal')
plt.sca(ax2)
plt.xlim([0, 200])
plt.ylim([0, 27])
plt.xlabel('Indoor rotational velocity (rad/s)', fontsize=12)
plt.ylabel('Outdoor distance traveled (m)', fontsize=12)
plt.savefig('IOrotv-dist.eps', transparent=True, bbox_inches='tight')
plt.savefig('IOrotv-dist.png', dpi=300, transparent=True, bbox_inches='tight')
plt.show()
Speeds = collected_data[collected_data['SamaraType']=='r']['IDspeedMean'].values
Rotvels = collected_data[collected_data['SamaraType']=='r']['IDrotMean'].values

plt.figure()
plt.scatter(Rotvels, Speeds, color='orange', edgecolor='black')
plt.xlabel('Rotational velocity (rad/s)')
plt.ylabel('Descent velocity (m/s)')

plt.figure()
plt.scatter(Rotvels*radii/1000, Speeds, color='orange', edgecolor='black')
plt.xlabel('Tip velocity (m/s)')
plt.ylabel('Descent velocity (m/s)')
plt.savefig('Figures/TipVelocity.png', dpi=300, transparent=True, bbox_inches=0)
from windrose import WindroseAxes
from matplotlib import patches

weatherData = pd.read_csv('DropStudy/OutdoorDropData/weatherData.csv')
weatherData['Times'] = weatherData['Timestamp'].apply(pd.Timestamp)
starttimes = ['2017-04-25T10','2017-04-25T10:40','2017-04-26T9:43','2017-04-27T9:42']
stoptimes = ['2017-04-25T10:20','2017-04-25T11:00','2017-04-26T10:43','2017-04-27T10:42']
df_list = [new_df[np.logical_and(new_df['Date']==pd.Timestamp('2017-04-24'), new_df['SamaraType']=='SM')], 
           new_df[np.logical_and(new_df['Date']==pd.Timestamp('2017-04-24'), new_df['SamaraType']=='NM')], 
           new_df[new_df['Date']==pd.Timestamp('2017-04-26')], 
           new_df[new_df['Date']==pd.Timestamp('2017-04-27')]]

fig, ax = plt.subplots(figsize=(3.25,3))
# fig2, ax2 = plt.subplots(nrows=1, ncols=3, figsize=(6.5,2.8), sharex=True, sharey=True)
fig8, ax8 = plt.subplots(nrows=2, ncols=2, figsize=(3.2, 3.2))
ax8 = [ax8[0,0], ax8[1,0], ax8[0, 1], ax8[1,1]]
fig3, ax3 = plt.subplots(figsize=(3.25,3))
fig4, ax4 = plt.subplots(nrows=1, ncols=3, figsize=(6.5,1.8), sharex=True, sharey=True)
fig5, ax5 = plt.subplots(figsize=(3.25, 3))
fig6, ax6 = plt.subplots(figsize=(3.25, 3))

# fig4, ax4 = plt.subplots(figsize=(1,1))
# fig5, ax5 = plt.subplots(figsize=(1,1))
# fig6, ax6 = plt.subplots(figsize=(1,1))
WindAxList = [WindroseAxes.from_ax() for k in range(4)]

NTimes = []
RTimes = []
STimes = []
NDists = []
RDists = []
SDists = []
# SDayDists = [[],[],[]]
# NDayDists = [[],[],[]]
# RDayDists = [[],[],[]]
STrialDists = [[], [], [], []]
NTrialDists = [[], [], [], []]
RTrialDists = [[], [], [], []]
STrialTimes = [[], [], [], []]
NTrialTimes = [[], [], [], []]
RTrialTimes = [[], [], [], []]

RWindage = []
NWindage = []
SWindage = []

for k, (start, stop, dropdata) in enumerate(zip(starttimes, stoptimes,df_list)):
    j = k
    start = pd.Timestamp(start)
    stop = pd.Timestamp(stop)
    temp=weatherData[np.logical_and(weatherData['Times']>start,weatherData['Times']<stop)]
    Directions = temp['Wind Vane'].values*np.pi/180
    Speeds = temp['Anemometer'].values*0.44704
    ax5.plot(Directions*180/np.pi)
    ax6.plot(Speeds)

### Windrose, day
#     if k == 0:
#         tempDirections = Directions.copy()
#         tempSpeeds = Speeds.copy()
#     elif k == 1:
#         tempDirections = tempDirections+Directions.copy()
#         tempSpeeds = tempSpeeds+Speeds.copy()
#         WindAxList[k-1].bar(tempDirections*180/np.pi, tempSpeeds, bins=[k for k in range(7)], blowto=True, normed=True, opening=0.95, edgecolor='white')
#         WindAxList[k-1].set_legend()
#     else:
### Windrose, time period
    WindAxList[k].bar(Directions*180/np.pi, Speeds, bins=[k for k in range(7)], blowto=True, normed=True, opening=0.95, edgecolor='white')
    
    WindAng = np.average(Directions, weights=Speeds)

    WindStd = np.sqrt(np.cov(Directions, aweights=Speeds))
#     WindAng = np.arctan2(np.mean(Speeds*np.sin(Directions)), np.mean(Speeds*np.cos(Directions)))
#     WindStd = np.std(np.arctan2(Speeds*np.sin(Directions), Speeds*np.cos(Directions))%(2*np.pi))
    WindSpeed = np.mean(Speeds)
    print(f'Wind trial {k}\nAngle: {WindAng*180/np.pi:3.1f} deg\nAngle std: {WindStd*180/np.pi:3.1f} deg\nSpeed: {WindSpeed:4.2f} m/s\n\n')
    
    # Wind angle gives clockwise angle with 0 = south. Measurement gives clockwise angle with 0=north.
    # Therefore, new angle = (180+WindAng-Ang)%360
    for row in dropdata.itertuples():
        if ~np.isnan(row.Velocity):
            angle = (np.pi+WindAng-row.Ang*np.pi/180)%(2*np.pi)
            if row.SamaraType == 'N':
                color = 0.9*np.array([[1,165/255,0]])
                marker = 'o'
                RTimes.append(row.Time*WindSpeed)
                RDists.append(row.Dist)
                RTrialDists[k].append(row.Dist*0.0254)
                RTrialTimes[k].append(row.Time)
            elif row.SamaraType == 'NM':
                color = np.array([[0,1,0]])
                marker = '^'
                NTimes.append(row.Time*WindSpeed)
                NDists.append(row.Dist)
                NTrialDists[k].append(row.Dist*0.0254)
                NTrialTimes[k].append(row.Time)
            else:
                color = np.array([[0.8,0.8,0]])
                marker = 'v'
                STimes.append(row.Time*WindSpeed)
                SDists.append(row.Dist)
                STrialDists[k].append(row.Dist*0.0254)
                STrialTimes[k].append(row.Time)
            ax.scatter(row.Time*WindSpeed, row.Dist*0.0254, marker=marker, edgecolors='k', s=16, c=color)
            ax3.scatter(row.Time, row.Dist*0.0254, marker=marker, edgecolors='k', s=16, c=color)
#             ax2[j].scatter(row.Dist*np.sin(angle)*0.0254, row.Dist*np.cos(angle)*0.0254, marker=marker, edgecolors='k', s=16, c=color)
            ax8[k].scatter(-row.Dist*np.sin(angle)*0.0254, row.Dist*np.cos(angle)*0.0254, marker=marker, edgecolors='k', s=16, c=color)
            wedge = patches.Wedge((0,0), 6*WindSpeed, 90-WindStd*180/np.pi, 90+WindStd*180/np.pi, ec=(0.4, 0.4, 0.4), fc=(0.8, 0.8, 0.8), zorder=-100)
#             ax2[j].add_patch(wedge)
            wedge = patches.Wedge((0,0), 6*WindSpeed, 90-WindStd*180/np.pi, 90+WindStd*180/np.pi, ec=(0.4, 0.4, 0.4), fc=(0.8, 0.8, 0.8), zorder=-100)
            ax8[k].add_patch(wedge)
            for rad in [5, 15, 25]:
                circle = patches.Circle((0,0), rad, fc=None, ec=(0.2, 0.2, 0.2), zorder=-99)
                circle.set_fill(False)
#                 ax2[j].add_patch(circle)
                circle = patches.Circle((0,0), rad, fc=None, ec=(0.2, 0.2, 0.2), zorder=-99)
                circle.set_fill(False)
                ax8[k].add_patch(circle)
    ax.set_aspect('equal')
#     ax2[j].set_aspect('equal')
    ax8[k].set_aspect('equal')
#     ax2[j].set_xlim([-25, 25])
#     ax2[j].set_ylim([-5, 25])

#     print(np.mean(temp['Anemometer'].values), np.mean(temp['Wind Vane'].values))


for j in range(3):
#     print(NDayDists[j])
    plt.sca(ax4[j])
#     ax4[j].hist(SDayDists[j], color=np.array([0.8,0.8,0]))
#     ax4[j].hist(NDayDists[j], color=np.array([0,1,0]))
#     ax4[j].hist(RDayDists[j], color=0.9*np.array([1,165/255,0]))
    ax4[j].hist([STrialDists[j], NTrialDists[j], RTrialDists[j]], range=((0, 30)), bins=12,
                color=[np.array([0.8,0.8,0]),np.array([0,1,0]),0.9*np.array([1,165/255,0])], 
#                edgecolor=[np.array([0.8,0.8,0]),np.array([0,1,0]),0.9*np.array([1,165/255,0])])
                edgecolor=None)
    plt.xlabel('Distance (m)')
plt.sca(ax4[0])
plt.ylabel('Count')
plt.savefig('Figures/DistHist.eps', transparent=True, bbox_inches='tight')
plt.savefig('Figures/DistHist.png', dpi=300, transparent=True, bbox_inches='tight')
    
TimeLists = [RTimes, NTimes, STimes]
DistLists = [RDists, NDists, SDists]
Colors = [0.9*np.array([1,165/255,0]), np.array([0,1,0]), np.array([0.8,0.8,0])]

means, errs = [], []
for (Times, Dists, color) in zip(TimeLists, DistLists, Colors):
    x = np.array(Times)
    x = x.reshape(len(x), 1)
    y = np.array(Dists)*0.0254
    m, res = np.linalg.lstsq(x, y, rcond=None)[0:2]
    ss = np.sum((y-np.mean(y))**2)
    print(m, 1-res/ss)
    ax.plot([0, 30], [0, m*30], color=color, zorder=-99)
    Coeffs = [D/T for (T,D) in zip(x, y)]
    means.append(np.mean(Coeffs))
    errs.append(np.std(Coeffs))
ax.plot([0, 30], [0, 30], '--', color=(0.4, 0.4, 0.4), zorder=-100)
ax.set_xlim((0, 30))
# print(means)
# print(errs)

fig7, ax7 = plt.subplots()
ax7.bar([0,1,2], means, yerr=errs, align='center', color=Colors, ecolor='black', capsize=10)
ax7.plot([-0.5, 2.5], [1, 1], '--', color=(0.4, 0.4, 0.4), zorder=-100)
ax7.set_xlim((-0.5, 2.5))
ax7.set_xticks([0,1,2])
ax7.set_xticklabels(['Natural Norway', 'Artificial Norway', 'Artificial Silver'])
plt.ylabel('Windage coefficient')
plt.savefig('Figures/WindageBarplot.eps', transparent=True, bbox_inches='tight')
plt.savefig('Figures/WindageBarplot.png', dpi=300, transparent=True, bbox_inches='tight')
    
plt.sca(ax)
plt.xlabel('Passive tracer distance', fontsize=12)
plt.ylabel('Distance traveled', fontsize=12)
plt.ylim([0, 25])
plt.savefig('Figures/windage.eps', transparent=True, bbox_inches='tight')
plt.savefig('Figures/windage.png', dpi=300, transparent=True, bbox_inches='tight')

# plt.sca(ax2[1])
# plt.xlabel('Crosswind landing', fontsize=12)
# plt.sca(ax2[0])
# plt.ylabel('Downwind landing', fontsize=12)
# plt.savefig('Figures/landing.eps', transparent=True, bbox_inches='tight')
# plt.savefig('Figures/landing.png', dpi=300, transparent=True, bbox_inches='tight')

plt.sca(ax3)
plt.xlabel('Descent time', fontsize=12)
plt.ylabel('Distance', fontsize=12)
plt.savefig('Figures/time-dist.eps', transparent=True, bbox_inches='tight')
plt.savefig('Figures/time-dist.png', dpi=300, transparent=True, bbox_inches='tight')

plt.sca(ax5)
plt.xlabel('Time', fontsize=12)
plt.ylabel('Wind angle (deg)', fontsize=12)
plt.savefig('Figures/wind-angle.eps', transparent=True, bbox_inches='tight')
plt.savefig('Figures/wind-angle.png', dpi=300, transparent=True, bbox_inches='tight')

plt.sca(ax6)
plt.xlabel('Time', fontsize=12)
plt.ylabel('Wind speed (m/s)', fontsize=12)
plt.savefig('Figures/wind-speed.eps', transparent=True, bbox_inches='tight')
plt.savefig('Figures/wind-speed.png', dpi=300, transparent=True, bbox_inches='tight')


for k in range(4):
    plt.sca(WindAxList[k])
    plt.gca().set_xticklabels([])
    plt.gca().set_yticklabels([])
    plt.savefig(f'Figures/windrose_day{k}.eps', transparent=True, bbox_inches='tight')
    plt.savefig(f'Figures/windrose_day{k}.png', dpi=300, transparent=True, bbox_inches='tight')

for k in range(4):
    plt.sca(ax8[k])
    if k in [0,2]:
        plt.xlim([-6.25, 6.25])
        plt.ylim([-1.25, 6.25])
    else:        
        plt.xlim([-25, 25])
        plt.ylim([-5, 25])
plt.savefig('Figures/landing-4panel.eps', transparent=True, bbox_inches='tight')
plt.savefig('Figures/landing-4panel.png', dpi=300, transparent=True, bbox_inches='tight')
for (stype, trialtimes) in zip(['Artificial silver', 'Artificial Norway', 'Natural Norway'], [STrialTimes, NTrialTimes, RTrialTimes]):
    for t in range(4):
        if len(trialtimes[t])>0:
            print(f'Trial {t} {stype} time: {np.mean(trialtimes[t]):4.2f}')
namelist = ['Artificial silver', 'Artificial Norway', 'Natural Norway']
colorlist = [np.array([[0.8,0.8,0]]), np.array([[0,1,0]]), 0.9*np.array([[1,165/255,0]])]
distlista = [STrialDists[:2], NTrialDists[:2], RTrialDists[:2]]
distlistb = [STrialDists[2:], NTrialDists[2:], RTrialDists[2:]]
for (stype, color, dista, distb) in zip(namelist, colorlist, distlista, distlistb):
    templist = []
    for t in range(2):
        dists = dista[t] + distb[t]
        print(f'Condition {t+1} {stype} distance: {np.mean(dists):4.2f} m')
        plt.figure()
        weights = np.ones((len(dists),))*1/len(dists)
        if t == 0:
            bins = 24
            xran = [0, 6]
        else:
            bins = 25
            xran = [0, 25]
        plt.hist(dists, weights=weights, bins=bins, range=xran, color=color)
        plt.xlim(xran)
        plt.ylim([0, 0.25])
        plt.xlabel('Distance (m)')
        plt.savefig(f'./Figures/hist-{stype}-group{t+1}.eps', transparent=True, bbox_inches=0)
        plt.savefig(f'./Figures/hist-{stype}-group{t+1}.png', dpi=300, transparent=True, bbox_inches=0)
ws = np.linspace(1.25, 3.25, 100)
for df in np.arange(1.5,5.5,0.5):
    plt.plot(ws, np.sqrt(0.61*df/ws), color=(0.7-0.1*df)*np.array([1,1,1]), zorder=-2)

wing_loadings = np.array([data_list[66+k][4] for k in range(30)])
Speeds = np.array(collected_data[collected_data['SamaraType']=='r']['IDspeedMean'].values)
# plt.scatter(wing_loadings, 1/Speeds, s=30, color='orange', edgecolor='black', zorder=0)
plt.errorbar(np.nanmean(wing_loadings), 1/np.nanmean(Speeds), xerr=np.nanstd(wing_loadings), yerr=np.nanstd(1/Speeds), color='orange', zorder=1)
plt.scatter(np.nanmean(wing_loadings), 1/np.nanmean(Speeds), s=60, color='orange', edgecolor='black', zorder=2)
Speeds2 = np.array(collected_data[collected_data['SamaraType']=='n']['IDspeedMean'].values)
plt.errorbar(3.183, 1/np.nanmean(Speeds2), xerr=.131, yerr=np.nanstd(1/Speeds2), color='green', zorder=1)
plt.scatter(3.183, 1/np.nanmean(Speeds2), s=60, marker='v', color='green', edgecolor='black', zorder=2)
plt.xlabel('Wing loading')
plt.ylabel('Descent time per meter (m/s)')
plt.savefig('Figures/LentinkFigure-var.eps', bbox_inches=0, transparent=True)
plt.savefig('Figures/LentinkFigure-var.png', dpi=300, bbox_inches=0, transparent=True)
ws = np.linspace(1.25, 3.25, 100)
for df in np.arange(1.5,5.5,0.5):
    plt.plot(ws, 1/np.sqrt(0.61*df/ws), color=(0.7-0.1*df)*np.array([1,1,1]), zorder=-2)

wing_loadings = np.array([data_list[66+k][4] for k in range(30)])
Speeds = np.array(collected_data[collected_data['SamaraType']=='r']['IDspeedMean'].values)
# plt.scatter(wing_loadings, Speeds, s=30, color='orange', edgecolor='black', zorder=0)
plt.errorbar(np.nanmean(wing_loadings), np.nanmean(Speeds), xerr=np.nanstd(wing_loadings), yerr=np.nanstd(1/Speeds), color='orange', zorder=1)
plt.scatter(np.nanmean(wing_loadings), np.nanmean(Speeds), s=60, color='orange', edgecolor='black', zorder=2)
Speeds2 = np.array(collected_data[collected_data['SamaraType']=='n']['IDspeedMean'].values)
plt.errorbar(3.183, np.nanmean(Speeds2), xerr=.131, yerr=np.nanstd(1/Speeds2), color='green', zorder=1)
plt.scatter(3.183, np.nanmean(Speeds2), s=60, marker='v', color='green', edgecolor='black', zorder=2)
plt.xlabel('Wing loading')
plt.ylabel('Descent speed (m/s)')
plt.savefig('Figures/LentinkFigure-invert-var.eps', bbox_inches=0, transparent=True)
plt.savefig('Figures/LentinkFigure-invert-var.png', dpi=300, bbox_inches=0, transparent=True)
import cv2
import skimage.morphology as skmorph

rfile = './ExampleVideos/r-g03-t03_20170719_153528_C001H001S0001.avi'
nfile = './ExampleVideos/n-g03-t03_20170719_153749_C001H001S0001.avi'

rdatfile = './DropStudy/VideoScripts/Real Samara Data/r-g03-t03-data.csv'
ndatfile = './DropStudy/VideoScripts/Norway Trial Data/n-g03-t03-data.csv'

rtimes = [315, 450, 590, 730, 905]
ntimes = [95, 195, 325, 445, 570]

def getimages(vid, times):
    frames = []
    for _ in range(times[0]):
        ok, frame = vid.read()
    frames.append(frame.copy())
    for time in np.diff(times):
        for _ in range(time):
            ok, frame = vid.read()
        frames.append(frame.copy())
    return frames

rvid = cv2.VideoCapture(rfile)
nvid = cv2.VideoCapture(nfile)

rframes = getimages(rvid, rtimes)
nframes = getimages(nvid, ntimes)

rvid.release()
nvid.release()
def plotimages(frames, bbox):
    fig, ax = plt.subplots(1, 5, figsize=(20, 10))
    for axi, frame in zip(ax.ravel(), frames):
        axi.imshow(frame[bbox[0]:bbox[2], bbox[1]:bbox[3]])
    return fig, ax

rbox = [50, 550, 450, 750]
nbox = [0, 850, 400, 1050]

rfig, rax = plotimages(rframes, rbox)
nfig, nax = plotimages(nframes, nbox)
def makecomposite(frames, bbox):
    composite = np.zeros((bbox[2]-bbox[0], bbox[3]-bbox[1]))
    for frame in frames:
        crop = frame[bbox[0]:bbox[2], bbox[1]:bbox[3], 0]
        composite += crop.astype(float)
    composite = composite/len(frames)
    for j in range(2):
        composite = composite-np.mean(composite)
        composite[composite<0] = 0
    return skmorph.opening(composite)

rcomp = makecomposite(rframes, rbox)
ncomp = makecomposite(nframes, nbox)

fig, ax = plt.subplots(1, 2)
ax[0].imshow(rcomp, cmap='gray')
ax[1].imshow(ncomp, cmap='gray')
ax[0].axis('off')
ax[1].axis('off')
import pandas as pd
import scipy.ndimage.filters as filt

def smoothpath(path, bbox, sigma = 5, cutoff=5):
    x = path['Column'].values
    y = path['Row'].values
    
    xsmooth = filt.gaussian_filter1d(x, sigma)    
    xnan = np.where(np.abs(x-xsmooth)<cutoff, x, np.nan)
    
    xsmoother = filt.gaussian_filter1d(xnan, 1)-bbox[1]
    ysmoother = filt.gaussian_filter1d(y, 1)-bbox[0]
    
    intfun = interp.interp1d(ysmoother[~np.isnan(xsmoother)], xsmoother[~np.isnan(xsmoother)], kind='linear')
    yout = np.linspace(ysmoother[0], ysmoother[-1], 1000)
    xout = intfun(yout)
    
    return xout, yout

rpath = pd.read_csv(rdatfile)
npath = pd.read_csv(ndatfile)

rx, ry = smoothpath(rpath, rbox)
nx, ny = smoothpath(npath, nbox)

fig, ax = plt.subplots(1, 2)
ax[0].imshow(rcomp, cmap='gray')
ax[0].plot(rx, ry, color='cyan')
ax[0].set_xlim((0, rbox[3]-rbox[1]))
ax[0].set_ylim((0, rbox[2]-rbox[0]))
ax[0].invert_yaxis()
ax[0].axis('off')

ax[1].imshow(ncomp, cmap='gray')
ax[1].plot(nx, ny, color='cyan')
ax[1].set_xlim((0, nbox[3]-nbox[1]))
ax[1].set_ylim((0, nbox[2]-nbox[0]))
ax[1].invert_yaxis()
ax[1].axis('off')
plt.savefig('Figures/Trajectory-plots.eps', transparent=True, bbox_inches=0)
plt.savefig('Figures/Trajectory-plots.png', dpi=300, transparent=True, bbox_inches=0)
# -*- coding: utf-8 -*-
"""
Created on Mon Mar 26 10:05:57 2018

@author: Nathaniel
"""

import pandas as pd
import numpy as np
from numpy import ma
import scipy.interpolate as interp
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
sns.set_style('ticks')
mpl.rc('mathtext', rm='serif', fontset='cm')
mpl.rc('font', family='serif', serif='Times New Roman')
labelfont = {'fontsize':9, 'fontname':'Times New Roman'}
tickfont = {'fontsize':8, 'fontname':'Times New Roman'}




df = pd.read_excel('Density Study/Samara Segment Mass.xlsx', sheet_name='Sheet1')
data_list = df.values.tolist()

#Initialize variables and arrays
mass_length = matrix = np.zeros((30,21))
mass = matrix = np.zeros((30,21))
length = matrix = np.zeros((30,21))
area = matrix = np.zeros((30,1))
center_mass = matrix = np.zeros((30,1))
center_percent = matrix = np.zeros((30,1))
total_data = matrix = np.zeros((30,200))
total_data2 = matrix = np.zeros((30,200))
total_data3 = matrix = np.zeros((30,200))
total_data4 = matrix = np.zeros((30,200))
total_data5 = matrix = np.zeros((2,200))
area2 = matrix = np.zeros((30,1))
colors = [mpl.cm.viridis(a) for a in np.linspace(0, 1, 30)]
mass2 = matrix = np.zeros((2,21))
length2 = matrix = np.zeros((2,21))
mass_length2 = matrix = np.zeros((2,21))
center_mass2 = matrix = np.zeros((2,1))
center_percent2 = matrix = np.zeros((2,1))
center_percent2 = matrix = np.zeros((2,1))

radii = np.zeros(30,)

#Plot mass fraction vs length fraction
for a in range(0,30):
    mass[a, 1:] = data_list[a+1][1:21]
    length[a, 1:] = data_list[a+33][1:21]
    mass_length[a] = [b*c for b,c in zip(mass[a],np.cumsum(length[a]))]
    center_mass[a] = np.nansum(mass_length[a])/np.nansum(mass[a])
    center_percent[a] = center_mass[a]/np.nansum(length[a])
    radii[a] = np.nansum(length[a, :])*(1-center_percent[a])
    X = ma.masked_invalid(length[a])
    Y = ma.masked_invalid(mass[a])
    Xperc = np.cumsum(X[~X.mask])/np.sum(X[~X.mask])
    Yperc = Y[~X.mask]/np.sum(Y[~X.mask])
    f = interp.UnivariateSpline(Xperc, Yperc, k=3, s=0)
    total_data[a] = f(np.linspace(0,.99,200))
    plt.figure(1)
    plt.title("Mass Fraction vs Length Fraction for Natural Samaras (fit in red)")
#    plt.plot(np.linspace(0, 0.99, 200), f(np.linspace(0, 0.99, 200)))
    plt.xlabel('Fraction of length', **labelfont)
    plt.ylabel('Mass fraction, $m$', **labelfont)
    plt.xticks([0, 0.25, 0.5, 0.75, 1.0], **tickfont)
    plt.yticks([0, 0.25, 0.5, 0.75, 1.0], **tickfont)
    plt.ylim([0, 1])
    plt.plot(np.linspace(0, 0.99, 200), f(np.linspace(0, 0.99, 200)), color=colors[a])

for a, c in enumerate(center_percent):
    plt.figure(2)
    plt.title("Center of Mass at Length Fraction for Natural Samaras (average in red)")
    plt.xlabel("Fraction of Length")
    plt.ylabel("Center of Mass")
    plt.xlim([0,1])
    plt.xticks([0, 0.25, 0.5, 0.75, 1.0], **tickfont)
    plt.yticks([0, 0.25, 0.5, 0.75, 1.0], **tickfont)
    plt.plot([c, c], [0, 1], color=colors[a])
    plt.plot([np.mean(center_percent),np.mean(center_percent)],[0,1], color='r')
  


#Plot mass vs length fraction
for a in range(0,30):
    mass[a, 1:] = data_list[a+1][1:21]
    length[a, 1:] = data_list[a+33][1:21]
    mass_length[a] = [b*c for b,c in zip(mass[a],np.cumsum(length[a]))]
    center_mass[a] = np.nansum(mass_length[a])/np.nansum(mass[a])
    center_percent[a] = center_mass[a]/np.nansum(length[a])
    X = ma.masked_invalid(length[a])
    Y = ma.masked_invalid(mass[a])
    g = interp.UnivariateSpline(np.cumsum(X[~X.mask])/np.sum(X[~X.mask]), (Y[~X.mask]), k=3, s=0)
    total_data2[a] = g(np.linspace(0,.99,200))
    plt.figure(3)
    plt.title("Mass vs Length Fraction for Natural Samaras")
    plt.plot(np.linspace(0, 0.99, 200), g(np.linspace(0, 0.99, 200)), color=colors[a])
    plt.xlabel('Fraction of length', **labelfont)
    plt.ylabel('Mass, $mg$', **labelfont)
    plt.xticks([0, 0.25, 0.5, 0.75, 1.0], **tickfont)


#Plot Density vs length fraction
for a in range(0,30):
    mass[a, 1:] = data_list[a+1][1:21]
    length[a, 1:] = data_list[a+33][1:21]
    area[a] = data_list[a+65][1]
    mass_length[a] = [b*c for b,c in zip(mass[a],np.cumsum(length[a]))]
    center_mass[a] = np.nansum(mass_length[a])/np.nansum(mass[a])
    center_percent[a] = center_mass[a]/np.nansum(length[a])
    X = ma.masked_invalid(length[a])
    Y = ma.masked_invalid(mass[a])
#    Z = ma.masked_invalid(area[a])
    h = interp.UnivariateSpline(np.cumsum(X[~X.mask])/np.sum(X[~X.mask]), (Y[~X.mask])/(area[a]), k=3, s=0)
    total_data3[a] = h(np.linspace(0,.99,200))
    plt.figure(4)
    plt.title("Area Density vs Length Fraction for Natural Samaras")
    plt.plot(np.linspace(0, 0.99, 200), h(np.linspace(0, 0.99, 200)), color=colors[a])
    plt.xlabel('Fraction of length', **labelfont)
    plt.ylabel('Area Density, $mg/mm^2$', **labelfont)
    plt.xticks([0, 0.25, 0.5, 0.75, 1.0], **tickfont)



#Plot linear density vs length fraction
for a in range(0,30):
    mass[a, 1:] = data_list[a+1][1:21]
    length[a, 1:] = data_list[a+33][1:21]
    mass_length[a] = [b*c for b,c in zip(mass[a],np.cumsum(length[a]))]
    center_mass[a] = np.nansum(mass_length[a])/np.nansum(mass[a])
    center_percent[a] = center_mass[a]/np.nansum(length[a])
    X = ma.masked_invalid(length[a])
    Y = ma.masked_invalid(mass[a])
    j = interp.UnivariateSpline(np.cumsum(X[~X.mask])/np.sum(X[~X.mask]), ((Y[~X.mask])/(X[~X.mask]))*np.sum(X[~X.mask])/np.sum(Y[~X.mask]), k=3, s=0)
    total_data4[a] = j(np.linspace(0,.99,200))
    plt.figure(5)
    plt.title("Linear Density vs Length Fraction for Natural Samaras (fit in red)")
    plt.plot(np.linspace(0, 0.99, 200), j(np.linspace(0, 0.99, 200)), color=colors[a])
    plt.xlabel('Fraction of length', **labelfont)
    plt.ylabel('Linear Density', **labelfont)
    plt.xticks([0, 0.25, 0.5, 0.75, 1.0], **tickfont)
    area2[a] = np.sum((((Y[~X.mask])/(X[~X.mask]))*np.sum(X[~X.mask])/np.sum(Y[~X.mask]))*((X[~X.mask])/np.sum(X[~X.mask])))



for a in range(0,2):
    mass2[a, 1:] = data_list[a+97][1:21]
    length2[a, 1:] = data_list[a+100][1:21]
    mass_length2[a] = [b*c for b,c in zip(mass2[a],np.cumsum(length2[a]))]
    center_mass2[a] = np.nansum(mass_length2[a])/np.nansum(mass2[a])
    center_percent2[a] = center_mass2[a]/np.nansum(length2[a])
    X2 = ma.masked_invalid(length2[a])
    Y2 = ma.masked_invalid(mass2[a])
    w = interp.UnivariateSpline(np.cumsum(X2[~X2.mask])/np.sum(X2[~X2.mask]), ((Y2[~X2.mask])/(X2[~X2.mask]))*np.sum(X2[~X2.mask])/np.sum(Y2[~X2.mask]), k=3, s=0)
    total_data5[a] = w(np.linspace(0,.99,200))
    plt.figure(6)
    plt.title("Linear Density vs Length Fraction for 3D Printed PLA Samaras (natural fit in red)")
    plt.plot(np.linspace(0, 0.99, 200), w(np.linspace(0, 0.99, 200)), color=colors[a])
    plt.xlabel('Fraction of length', **labelfont)
    plt.ylabel('Linear Density', **labelfont)
    plt.xticks([0, 0.25, 0.5, 0.75, 1.0], **tickfont)
    area2[a] = np.sum((((Y[~X.mask])/(X[~X.mask]))*np.sum(X[~X.mask])/np.sum(Y[~X.mask]))*((X[~X.mask])/np.sum(X[~X.mask])))
        
        

        
# Find best fit line for mass fraction vs length fraction
#area = np.trapz(total_data[5],dx=.1)
xdata = np.linspace(0,.99,200)
mean_data = [np.mean(total_data[:, k]) for k in range(len(xdata))]
mean_data2 = [np.mean(total_data4[:, k]) for k in range(len(xdata))]
std_data = [np.std(total_data[:, k]) for k in range(len(xdata))]
std_data2 = [np.std(total_data4[:, k]) for k in range(len(xdata))]
#xdata = np.vstack([xdata[1:199:1] for k in range(30)]).ravel()
xdata = np.vstack([xdata for k in range(30)]).ravel()
#remove 0 entry from array
#total_data = total_data[:,1:199:1].ravel()
total_data = total_data.ravel()
total_data4 = total_data4.ravel()

#fit data
z = np.polyfit(xdata,total_data,9)
z2 = np.polyfit(xdata,total_data4,9)

xdata = np.linspace(0,.99,200)
y = np.poly1d(z)
y2 = np.poly1d(z2)
plt.figure(1)
plt.plot(xdata, y(xdata), color='r')
plt.figure(5)
plt.plot(xdata, y2(xdata),color='r')
plt.figure(6)
plt.plot(xdata, y2(xdata),color='r')
#plot mean data
#plt.plot(np.linspace(0,.99,200), mean_data, color='k')



#Plot fit by itself
#plt.figure(5)
#plt.plot(xdata,y(xdata), color='r')
#plt.xlabel('Fraction of length', **labelfont)
#plt.ylabel('Mass fraction, $m$', **labelfont)

plt.show()
#area = np.trapz(total_data[3],dx=.1)

plt.figure(7)
ax = plt.gca()
ax.errorbar(xdata, mean_data, yerr=std_data)
ax.plot(xdata,y(xdata), color='r')
plt.title("Standard Deviation for Mass Fraction Fit")

plt.figure(8)
ax = plt.gca()
ax.errorbar(xdata, mean_data2, yerr=std_data2)
ax.plot(xdata,y2(xdata), color='r')
plt.title("Standard Deviation for Linear Density Fit")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1995745812283967, 'fraction_numerical': 0.04974740760436054, 'mean_word_length': 4.14584758516897, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '124312', 'license': 'MIT', 'n_tokens_mistral': 16400, 'n_tokens_neox': 15154, 'n_words': 2416, 'path': 'AllDataViz.ipynb', 'repo_name': 'RossDynamics/team-samara', 'size': '37610'}