function [ im_main ] = createmainimage(Xus,r,theta,num_pixel)

%%%this code block is used to create the images %%%%%
%%%later part of this code block is used to save images %%%%%%%

%%theta=(2*pi)/180; %in radian
%%r=2;
%%Xus=10;
angle=-(pi/10):theta:(pi/10);  %%outf this angle no images are taken 
[alpha,beta]=size(angle);
%num_pixel=Xus/0.1;

im=zeros(Xus/0.1,Xus/0.1,beta);
for alpha=1:beta
    notice=waitbar(0,'loading visualization at an angle, please wait...');
    % Define the input grid (in 3D) %%that sutta code
    figure;
    
    [x3, y3, z3] = meshgrid(linspace(-10,10));
    % Compute the implicitly defined function x^2 + (y-6)^2 + z^3 - 2^2 = 0
    f1 = x3.^2 + (y3-6).^2 + z3.^2 - 2^2;                    %%why is this taking an elipse 

    % This surface plane, which can also be expressed as
    f2 = z3-tan(angle(alpha))*y3 - 0*x3.^3 ;
    % Also compute plane in the 'traditional' way.
    [x2, y2] = meshgrid(linspace(-10,10));
    z2 = tan(angle(alpha))*y2 - 0*x2;
    % Visualize the two surfaces.
    patch(isosurface(x3, y3, z3, f1, 0), 'FaceColor', [0.5 1.0 0.5], 'EdgeColor', 'none');
    patch(isosurface(x3, y3, z3, f2, 0), 'FaceColor', [1.0 0.0 0.0], 'EdgeColor', 'none');
    view(3); camlight; axis vis3d;
    set(gca,'xlim',[-5 5 ], 'ylim',[1 11], 'zlim',[-5 5]);
    xlabel('X');ylabel('Y');zlabel('Z');
    
    %title('visualization at',angle(alpha),'radians');
    
    
    % Find the difference field.
    f3 = f1 - f2;
    % Interpolate the difference field on the explicitly defined surface.
    f3s = interp3(x3, y3, z3, f3, x2, y2, z2);

    % Find the contour where the difference (on the surface) is zero.
    C = contours(x2, y2, f3s, [0 0]);
    % Extract the x- and y-locations from the contour matrix C.
    xL = C(1, 2:end);
    yL = C(2, 2:end);
    % Interpolate on the first surface to find z-locations for the intersection
    % line.
    
    if ~isempty(xL) & ~isempty(yL) 
        zL = interp2(x2, y2, z2, xL, yL);
        % Visualize the line.
        line(xL,yL,zL,'Color','k','LineWidth',3);

        %xL_1=ceil(xL/0.2)+50; %conversion to pixels
        %yL_1=ceil(yL/0.2)+50; %conversion to pixels   %%commented since errors

        max_xL=max(xL);
        min_xL=min(xL);
        max_yL=max(yL);
        min_yL=min(yL);

        
        %%direct substuition into eclipse formula
        a=(max_xL-min_xL)/2; %average to reduce errors
        b=(max_yL-min_yL)/2; 

        a=a/0.1;        %convertion to pixels
        b=b/0.1;

        H=0;Y=0; %center of the eclipse
        for i=1:num_pixel
            for j=1:num_pixel
                    if ((i-num_pixel/2)/a)^2+((j-num_pixel/2)/b)^2 <= 1 %no need to worry about the world cordinates for now
                        im(i,j,alpha)=255;
                    end
            end
        end
    end
    
%figure, 
%imshow(im);
    close(notice)
end


%%%%%%image save%%%%%%%%

%%% a total of 61 image slices will be saved <as 0 is included as well>%%%%%%%%
%%% need to add a total of 21 empty images both at the front and back %%%%%%%%% 

im_main=zeros(Xus/0.1,Xus/0.1,61);
for i=[22:40]
    im_main(:,:,i)=im(:,:,i-21);
end

%%run this loop if you wish to save the images
 %%image_index=1;
 %%filename = 'p%d.bmp';
 %%filename = sprintf(filename,image_index);
 %%image_index = image_index + 1;
 %%imwrite(im,filename);
    
im_main=uint8(im_main);
end