% Measurement Of A Gravity Anomaly Due To A Rising Plume Head

density = 5.0;  % g/cm^3
G = 6.67*10^-8; % dyne cm^2/g^2
k = 0.005;
Fe = 50000;
theta = 60;
jMax = 480;     % variable for calculating correct x-coordinate
f = 0;          % movie frame variable

for i = 1:5:450 
    for j = 1:960 % (-480 to +480)
        
        if j <= 480
            jj = -1 * jMax + j; % negative x-values (-480-0)
        else jj = j - jMax;     % positive x-values (1-480)
        end
        
        %plume head
        imageSizeY = 640; % x-window dimension
        imageSizeX = 480; % y-window dimension
        [columns, rows] = meshgrid(1:imageSizeX, 1:imageSizeY); %generate 2D rectangular grid
        centerX = 240;      % x-position
        centerY = 501-i;    % y-position
        radius = 50;        %radius of plume head
        circle = (rows - centerY).^2 + (columns - centerX).^2 <= radius.^2;   % calculate size of plume head

        gravity = ((G*4*pi*(radius^3)*density)/3)*(centerY/((jj^2 + centerY^2)^(3/2))); %calculate gravity anomaly
        
        % vertical magnetic anomaly
        Za = ((((4/3)*pi*radius^3*k*Fe)*sind(theta))/(jj^2+centerY^2)^(3/2)) * ((3*centerY^2)/(jj^2 + ...
            centerY^2)) - (((3*jj*centerY)/(jj^2 + centerY^2))*(1/tand(theta))-1);
        
%         % horizontal magnetic anomaly
        Ha = ((((4/3)*pi*radius^3*k*Fe)*cosd(theta))/(jj^2 + centerY^2)^(3/2)) * ((((3*jj^2)/(jj^2 + ...
            centerY^2)) - 1) - (((3*jj*centerY)/(jj^2 + centerY^2))*tand(theta)));
        
        gAnomaly(j,1) = gravity;    % store gravit anomaly values
        
        Zmag(j,1) = Za;             % store vertical magnetic anomlay values
        
        Hmag(j,1) = Ha;             % store horizontal magnetic anomaly values
        
        totalMag = Zmag*sind(theta) + Hmag*cosd(theta);     % total magnetic anomaly values
        
        xValue(j,1) = jj;           % store x-coordinate
        
    end
    
        f = f + 1;  % movie frame counter
        
%         figure(1)
%         subplot(1, 2, 1)    
%         image(circle) ;                 % plot rising plume head
%         colormap([0.8 0.4 0; 1 0 0]);
%         title('Rising Plume Head');
%         xlabel('Distance (m)')
%         ylabel('Depth (m)')
%         axis([0 480 0 640])
        
%         subplot(1, 2, 2)
%         plot(xValue, gAnomaly, '-');    % plot gravity anomaly
%         title('Gravity Anomaly')
%         xlabel('Distance (m)')
%         ylabel('Gravity Anomaly (mGal)')
%         axis([-480 480 0 0.00006])
        
        figure(2)
        plot(xValue, Zmag, '-', xValue, Hmag, 'r--', xValue, totalMag, 'g-');       %plot vertical & horizontal magnetic anomaly
        title('Magnetic Anomaly')
        xlabel('Distance (m)')
        ylabel('Magnetic Anomlay (nT)')
        legend('Horizontal Component', 'Vertical Component', 'Total Anomaly')
        
        saveas(gcf, ['magnetic_',sprintf('%03d',f),'.png'])
        pause(0.01)
end