Entering my senior year at school, I am looking ahead to my capstone course (winter 2021) and research opportunities at W&L. My goal is to develop a machine learning algorithm that is capable of accepting previously collected borehole data to predict the distribution and volume of oil and gas pockets within in an active oil field. Although I have had some exposure to developing machine learning algorithms in prior coursework, I am currently developing a financial forecasting model in an effort to become more comfortable with supervised learning programs. Because of the frequently updated data published by stock exchanges, I hope to be able to continuously train and test my model against current stock trends.

What are RR Lyrae Stars?

In the fall of 2020, I joined a research team studying RR Lyrae type variable stars. This class of star generally exhibits a periodic behavior that allows us to estimate luminosity, and therefore absolute magnitude according to the equation below:

$M$₁ - M₂ = -2.5 × log(^L₁ ∕ _L2)

where M₁ is the absolute magnitude of star 1, M₂ is the absolute magnitude of star 2, L₁ is the luminosity of star 1, and L₂ is the luminosity of star 2. If we examine the equation, we know that in order to calculate the absolute magnitude of star 1, we need to know its own lumonisty (easily found for RR Lyrae stars) and a second, reference stars' absolute magnitude and luminosity. So, by measuring the luminosity of a stellar object (a comparatively simple task) we are able to generate an estimate of the star’s absolute magnitude. We can then combine the absolute magnitude with our observed apparent magnitude to estimate the distance to the star as described by:

$d\; =\; 10(m\; -\; M\; -\; A\; +\; 5)\; ∕\; 5$

where m is the apparent magnitude, M is the absolute magnitude, d is the distance to the star, and A is the extinction factor. For our research, we performed the absolute magnitude and distance calculation several times with different filters to ensure more accurate results. We explain this process in more detail here.

Thus, because of a simple periodic behavior inherent in RR Lyrae stars, we are able to produce an accurate distance measurement! This unique property has landed RR Lyrae stars the designation of standard candles-- objects that serve a distance reference for the surrounding stellar objects. Although perhaps relatively abstract, this places special significance on this class of stars. If we are able to locate an RR Lyrae star in a galaxy, galactic cloud, or any other celestial body, we can derive a distance measurement to that location and use it as a jumping point for increasingly farther distance calculations. With a network of such standard candles in place, we can begin charting a huge number of systems that currently have little or no recorded information. Of course, this process becomes much more challenging when applying theory to reality, but it still remains a particularly accurate and comparatively easy method for estimating distances.

Our Role with OZ UMa

Although the project is based out of Australia, we were able to conduct our work remotely in coordination with our Washington and Lee-based research advisor, David Sukow.

As members of the larger research group, our role was to select a star and carry out the aforementioned process for estimating distance. My team, consisting of two W&L classmates, selected OZ UMa has our target. OZ UMa is an RRab star (a sub-class RR Lyrae stars) situated in Ursa Major and has very little information available from its limited previous studies. Using publicly available telescope time through Las Cumbres Observatory Network, we tracked OZ UMa over a 14 night period to gather information on the stars period fluctuation. Our efforts have culminated in a number of new or updated measurements of OZ UMa’s stellar properties. Although currently awaiting publication to the Astronomy Theory, Observations, and Methods Journal, you can read more about OZ UMa and our research through this pdf.

432 Park Ave. is one of New York's City most recent additions to the city skyline. Completed in 2015, the estimated cost of the project sits between a staggering $1.25 and $3.12 billion. Although construction costs are partially attributable to its location on Billionaire’s Row, the building’s innovative design certainly influenced the astronomical cost. Measuring at a 15:1 height-to-width ratio, 432 Park Ave. is currently one of the most slender structures in the world placing it in a new class of buildings described as “Pencil Towers.”

Much like balancing a pencil vertically on its eraser, fashioning such a building requires an extreme amount of precision. But, for a pencil, you need not even blow on it to make it fall on its side-- so how does a building keep from collapsing when subjected to hurricane-force winds, earthquakes, the commotion of underground railroads, and its own internal reverberations? The simulation to the right (provided by Strand7) provides a good visualization of extreme movement. As it turns out, the solution is not to stop the building from moving, but rather to control it. For 432 Park Ave., this required innovative design procedures like “gap floors” that dilute the effects of wind.

This area of study is of extreme practical importance when designing structures and is the focus of our investigation into 432 Park Ave. Specifically, we study how the building is constructed to ensure that it is not forced to oscillate in its resonant frequency which can produce catastrophic side-effects.

Natural phenomenon like earthquakes and hurricane force winds are a common and potentially dangerous influence on a building’s movement. Because these events generally oscillate at very low frequencies, when modelling their potential effects, we can limit our scope of study to the first ten natural frequencies of a system. In our case, if we know the first 10 natural frequencies of 432 Park Ave., we can compare them with historical data and machine generated simulations of natural phenomenon. This allows engineers to ensure that the structure will not have a natural frequency that is equal to a frequency produced by the environment. The GIF to the left is an example of a bridge (the Tacoma Narrows Bridge) that was forced to resonate in its natural frequency by strong gusts of wind and later collapsed.

To model this system, we utilized a degree-of-freedom model where each floor added a single degree-of-freedom (see accompanying illustration to the right). Thus, with 96 total floors, including basement and “gap” floors, we have 96 degrees-of-freedom and, in turn, our model produces 96 coupled differential equations that describe the buildings motion. This allows us to consider the buildings support beams as springs, where the equation for the spring constant of a single support beam is given by:

$k\; =\; 48\; \times \; (E\; \times \; I)\; ∕\; L3$

Casting our degree-of-freedom equations as a matrix allows us to calculate the modes and mode shapes of the the structure by solving for the eigenvectors and eigenvalues. When solving this system of equations, the eigenvalues correspond to the building's natural frequencies, while the eigenvectors describe the motion (or mode shapes) of the building for a specific eigenvalue. Each eigenvalue will have a corresponding eigenvector to describe its motion. As you might imagine, this is not a particularly fun equation to write out by hand, so we turned to MatLab to perform this complex computation.

To to do this, we utilized and expanded on previous MatLab code to compute our values. The drop down below contains the original code and the modifications we used to generate the eigenvectors and eigenvalues of our system. It also contains the script for the corresponding simulation rendered in MatLab.

%% script for running N story building sim.  Two examples are included.   
% 1. No initial conditions specified. In this case, modes are equally mixed with all c_n's set equal and phases all set to 0. 
% 2. Initial conditions specified.  Note that you could use the eigIC.m script to solve for C and Phi first. 
% Last modified; April 2020
 
%% Example with no Initial conditions specified.  
%specifying single values for m and k, NstoryBulidingSim generates the appropriate M and K matrices. 
%Other option is to fully specify M and K as matrix and pass those the main workhorse function NStoryBuildingSim(...) 

% %% 2. Example specifying initial conditions 
% Nfloors = 2; 
% m = 2e5; 
% k = 4e5; 
% Phi = [pi/2, pi/2]; 
% C = 0.2*[0.9772; -2.5583]; 
% [Omega, V] = NStoryBuildingSim(Nfloors, m, k, C, Phi);

%Our Values
m = 487350; %kg 
k = 55600000000; %N/m 
Nfloors = 96; 
close all; 
[Omega, V] = NStoryBuildingSim(Nfloors, m, k); 


            -------------------------------------------------------------------------------------------------------------------------------
    
    

% [Omega, V, simtimes, Xout, f, Fspec] = NStoryBuildingSim(Nfloors, m, k c, phi)
% Simulates building modes with NFLOORS.  
% The mass of each floor is M, the total support/sping constant for each floor is K.
% Each floor and set of support columns is assumed identical
function [Omega, V, simtimes, Xout, f, Fspec] = NStoryBuildingSim(Nfloors,m,k, c, phi)

% Initial conditions can be input with C and PHI (Nfloors x 1). 
% Plots are generated showing the mode shapes swaying in time, the displacement of each floor.  When the simulation finishes the FFT is
% displayed. Outputs allow further processing and plotting, if desired. By default, the total simulation time is equal to 3 periods of the
% largest natural period (smallest wo)
%
%
%
% ---- INPUTS ----
% Nfloors - number of floors in the building. Default = 3.
% c - mode mixing constants (Nfloors x 1).  Default is equal mode mixing:
%       ones(1,Nfloors)/Nfloors, e.g. 1/3*[1 1 1] 
% phi - initial phase angles. Default: zeros.
%
% 
%
% ---- OUTPUTS-----
%
% Omega  - vector of natural frequencies (rad/s)
% V -  Modal matrix. Each column corresponds to a (normalized) mode shape
%
% simtimes - vector of times at which positions are computed
% Xout - X position of each floor vs. time. Data is rows.
%
% f - vector of frequencies at which frequency spectra FSPEC is computed.
% Fspec - frequency spectra. Data in rows.
%
% -------- EXAMPLE -----
% Nfloors = 4;
% c = [1, 0.5, 0.25, 0.75];
% NStoryBuildingSim(Nfloors, c);
%
% See also eig, diag
%
% Original Author: Jon Erickson, W&L, ericksonj@wlu.edu 
% Last updated: April 2020 

%% JE Mar 23 2010
clrwheel2 = cell2mat({'b'; 'r'; 'c';'g'; 'm'; 'k'});

Ncycles =3;
%% choose a plot style
plotoverlay = false;

%% choose initial condition coeffs;
if nargin<1 | isempty(Nfloors)
    Nfloors = 3;
end

if nargin<4 |isempty(c)
    c = ones(1,Nfloors)/Nfloors;
end

if nargin<5 | isempty(phi)
    phi = zeros(1,Nfloors);
end

%% define matrices and solve for eigvects and vals
%mass matrix
%assumes each floor has same mass
M = m*eye(Nfloors);

%floors 86 and above are heavier by factor of 1.89
M(86:end, :) = M(86:end, :)*1.89;
% M2 = m2*eye(Nfloors);

%stiffness matrix
%assumes all columns are the same

% if Nfloors < 86
K = diag(2*k*ones(1,Nfloors))+diag(-k*ones(Nfloors-1,1),1)+diag(-k*ones(Nfloors-1,1),-1);
K(Nfloors, Nfloors) = k;
% else
    
K(86:end, :) = K(86:end, :)*0.78;
    
% K2 = diag(2*k2*ones(1,Nfloors))+diag(-k2*ones(Nfloors-1,1),1)+diag(-k2*ones(Nfloors-1,1),-1);
% K2(Nfloors, Nfloors) = k2;
% K = [k1+k2 -k2 0 ;
%     -k2 k2+k3 -k3;
%     0 -k3 k3];

%dynamic matrix
 D =inv(M)*K;
% if Nfloors < 86
%     D =inv(M)*K;
% else Nfloors > 86
%     D2 = inv(M2)*K2;
% end

[Vall,LambdaAll] = eig(D); % V is eigenvectors in columns. D is diag matrix of eigenvals.
OmegaAll = sqrt(diag(LambdaAll )); %vector of natural frequencies (rad/s).

% grab the 10  (nmodes) loest frequencies of oscillation
Nmodes = 10;
[OmLow, IIom] = sort(OmegaAll, 'ascend');
 Omega = OmLow(1:Nmodes);
V = Vall(:, IIom);

%normalize eigenvectors
for j = 1:Nmodes
    V(:,j) = V(:,j)/norm(V(:,j));
end

%sort the modes according to frequency:
[Omega,I] = sort(Omega, 'ascend');
V = V(:,I);

whos V Omega c
%computer natural periods:
To = 2*pi./Omega;
DT = min(2*pi./Omega)/10;
%% plot eigenvectors/shapes.
modefig = figure(1);

% modefignum = modefig.Number;
% figure(modefig);
    set(modefig, 'units','normalized', 'position', [ 0.0469    0.2046    0.5917    0.6889]);

% set(gcf, 'position', fullfig());
xshift = 1;

% close all;
iter = 1;
Ncycles*max(2*pi./Omega);
simtimes = 0:DT:Ncycles*max(2*pi./Omega);  %JE 09 Mar 2018: correcte dtime interval
Xout = zeros(Nfloors,length(simtimes));

for t = simtimes
    %compute total response at time t
    %     legend(num2str(To), 'location', 'northeast')
    X = zeros(Nfloors,1);  %total response. We'll sum to find it
    whos X V omega c
    for j = 1:Nmodes
        X = X + c(j)*V(:,j)*cos(Omega(j)*t + phi(j));
    end
    Xout(:,iter) = X;
    %then plot
    offset = zeros(1,Nmodes);
   
    figure(modefig.Number);
%     set(modefig, 'units','normalized', 'position', [ 0.6469    0.5046    0.2917    0.3889]);
    for j = 1:Nmodes+1
        if j==1
            %             offset(j) = max(V(:,j));
            cla;
        else
            %             offset(j) = max(V(:,j))+max(V(:,j-1));
            hold on;
            grid on;
        end
        
        %% Plot either overlaid or separate modes
        
        if j<=Nmodes
            clrJ = 1+mod(j,length(clrwheel2));
            % this plot command overlays the modes as they are plotted.
            if plotoverlay
                plot([0;V(:,j)]*cos(Omega(j)*t+phi(j)), 0:Nmodes, [clrwheel2(clrJ),'o:'], 'markerfacecolor', clrwheel2(j));
            else
                plot([0;V(:,j)]*cos(Omega(j)*t+phi(j))+(j-1)*xshift, 0:Nfloors, [clrwheel2(clrJ),'o:'], 'markerfacecolor', clrwheel2(clrJ));
                %                 plot([0;V(:,j)]*cos(Omega(j)*t)+xshift*(j-1), 0:Nmodes, [clrwheel2(j),'o:'], 'markerfacecolor', clrwheel2(j));
            end
        else
            if plotoverlay, plot([0;X], 0:Nfloors,'ko-', 'markerfacecolor', 'k'); else
                plot([0;X]+xshift*(j-1), 0:Nfloors,'ko-', 'markerfacecolor', 'k');end
        end
    end
    
    xlabelcell = {};
    for nn = 1:Nfloors
        xlabelcell{nn} = num2str(nn);
    end
    xlabelcell{Nfloors+1} = 'Total';
    set(gca, 'xticklabel', xlabelcell, 'xtick', xshift*[0:Nmodes]);
    %     set(gca, 'xticklabel', {'1'; '2'; '3'; 'Total'}, 'xtick', xshift*[0:Nmodes]);
    
    if ~plotoverlay
        xlim(1.05*[-xshift, xshift*(Nmodes+1)]);
    else
        xlim(1.05*xshift*[-1 1]);
    end
    ylim([0 Nmodes+1]);
    title(sprintf('Time = %2.3f s', t));
    pause(0.33);
    if iter==1; pause; end
    
    %% plot the position of each mass over time.
    figure(2); 
    set(gcf, 'units', 'normalized', 'position', [ 0.6443    0.03    0.2917    0.85]);
    for kk = 1:Nmodes %top floor is actually last mode.
        if kk ==1, xlabel('Time (s)'); end
        subplot(Nmodes,1, Nmodes-kk+1); hold on; xlim([0 simtimes(end)]);
        plot(simtimes(iter), Xout(kk,iter), 'ko-', 'markersize', 4);
        ylabel(sprintf('x_%d', kk));
        grid on;
        xlabel('Time (s)')
        
    end
    
    iter = iter+1;
    
end
figure(modefig.Number);
wo = 2*pi./To;
legend(num2str(To), 'location', 'northeast')

%% plot the FFT for each floor
L = length(Xout);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(Xout',NFFT)/L;
Fs = 1/DT;
f = Fs/2*linspace(0,1,NFFT/2+1);

figure(3);
set(gcf, 'units', 'normalized', 'position', [0.3432    0.0417    0.2917    0.3889]);
% Plot single-sided amplitude spectrum.

for ii = 1:Nmodes
    Fspec(:,ii) = 2*abs(Y(1:NFFT/2+1,ii));
    subplot(Nmodes,1,Nmodes-ii+1); %plot top floor on top. bottom floor on bottom
    plot(f,Fspec(:,ii));
    axis tight
    % Fspec = abs(Y);
    title(sprintf('FFT of x%d(t)',ii))
    xlabel('Frequency (Hz)')
    ylabel('|FFT(f)|')
end
                                    

This process can only be described as a first order approximation at its very best. To keep this project from getting out of hand, we do not consider particular construction designs (namely outriggers) or the effects of tuned mass dampers (TMDs). Moreover, we approximate construction design and materials to the best of our knowledge. To read more about 432 Park Ave. and degree-of-freedom systems visit the links below:

• New York Times Article
• Wikipedia Page
• The Architect's Newspaper
• YouTube video of TMD in action
• Degree-of-freedom PDF
• Video explanation of degree-of-freem system

I created this ballistic motion calculator to serve as the final project for the Washington and Lee Physics Department Modelling and Simulations in Physics course. The intended (and successful!) outcome of the project was to merge my physics and coding knowledge to develop a program capable of calculating and displaying the trajectory of a projectile under user-provided conditions.

To provide a sense for the breadth of the project, listed below are some of the conditions the calculator is able to account for:

1. Elevation
2. Temperature
3. Exit velocity/muzzle velocity
4. Air density
5. Projectile mass
6. Projectile cross area
7. Projectile diameter
8. Projectile ballistic coefficient
9. Target distance
10. Wind speed and direction

These inputs were then supplied to my python program that performed the appropriate calculations to describe the projectiles motion. Drawing upon the Euler-Lagrange equations, my program is designed to update the trajectory of the projectile in discrete, small time steps. The picture below is a stillshot of the animation generated by the program as it follows the projectiles path in times steps of dt = 0.05. At the sacrifice of computation time, we can use increasingly smaller time steps to produce increasingly more accurate outputs. However, in this case, A time step of 0.05 is sufficient.

My first challenge was to acquire the relevant inputs. To do this, I used the tkinter GUI interface to prompt the user for original inputs. However, I was careful to design the input request to ensure that potential future updates would allow for these inputs to be provided by external sources. For example, a modified Raspberry Pi would capable of providing the majority of these inputs for real-time, accurate data.

Next, I had to determine the appropiate equations to account for variables such as drag, wind, elevation, etc. Because of the difficulty associated with accurately accounting for drag, I opted to request the user for the projectile's ballistic coefficient. Generally, this value is known or even listed by the manufacturer and is much better suited for performing calculations on a wide variety of projectiles that can behave quite uniquely. To consider the effects of wind, I merely needed to decompose the vector input into constiuent components (x,y, and z direction) and include these components in the equation to update velocity. This means that our equations for velocity will use the relative velocity of the projectile.

Once the variable equations are constructed, the most important step can be included; the main equations to update motion! As previously mentioned, I draw upon the Euler-Lagrange equations to create a system of second order differential equations that will update the acceleration, and thus the velocity and position, at a time step of 0.05 (units depend on the users' inputs).

Finally, we can animate the motion of our projectile! To do this, I utilize the MatPlotLib Python package. This allowed me to create a 3-D environment that tracks and plots the motion of the projectile with the appropiate forces influencing its motion.

A social network is a group of people, each of whom is connected with a subset of the others. Such a network is represented using vertices (to represent each person) that are joined by edges (representing some sort of connection).

These social networks have been a topic of interest in the social sciences for over fifty years. Models have been developed to predict the spread of disease through a group of people, the relative connectivity of random individuals, and a variety of other important aspects of social networks.

Our project was aimed at building a social network based on scientific collaboration at an elite liberal arts college. Specifically, we focused on published researchers in the sciences at Washington and Lee University. Two of these researchers are considered connected if they have ever authored a paper together. This parameter is derived from the canonical paper from The Santa Fe Institute, "The Structure of Scientific Collaboration Networks." We restricted our data to Washington and Lee’s “hard sciences” departments: biology, chemistry, geology, physics, and biochemistry. This is in response to the tendency of social science professors to abstain from collaboration.

We subsequently compare these findings to Newman’s conclusions, a well-known application of the theory regarding the average degree of connectivity between individuals. Our result establishes a connection between a small, elite liberal arts college and the research tendencies of its faculty members. Future research could be aimed at extending this research to other similar-sized institutions for comparison and analysis.

For a summary of our results, please follow this link.