function StarfieldSimulation
% Source code that produces a GIF.
% Inspired by the "Starfield Simulation" Screensaver for Windows.
%
% 2017-04-20 Jahobr (update 2019-02-07 Jahobr)

[pathstr,fname] = fileparts(which(mfilename)); % save files under the same name and at file location
 
uniqueDepth = 200;
nFrames = 240;
renderDist = 200; % should be smaller or equal then uniqueDepth
StarDensity = 0.0005;
imageMax = [400, 500]; % pixel (height width)
subsampPix = 10; % subsampling per pixel in x and y
starRenderLimInPixel = 0.4; % diameter at least a x% of a pixel % if too small it flickers

for currentCase = 1:2
    rng(1); % int rand generator to make GIF repeatable; If you dont like the GIF try a differnt value
    speed = uniqueDepth/nFrames;
    nStars = StarDensity*(2*renderDist * 2*renderDist * uniqueDepth);
    
    x0 = rand(1,nStars)*renderDist*2-renderDist;
    y0 = rand(1,nStars)*renderDist*2-renderDist;
    sizze = 60+rand(1,nStars)*60;
    sizze = [sizze sizze]; % copy space behind each other to loop the gif
    
    depth = rand(1,nStars)*uniqueDepth;
    depth = sort(depth); % makes later rendering more efficient
    depth = [depth depth+uniqueDepth]; % copy space behind each other to loop the gif
    
    switch currentCase
        case 1 % classic
            saveName = fname;
            x_y =     [x0 x0;y0 y0]; % copy space behind each other to loop the gif
            thetaList = zeros(1,nFrames); % rotation
        case 2 % rot 90° per gif loop
            saveName = [fname '_rot'];
            x_y =     [x0 y0;y0 -x0]; % copy space behind each other to loop the gif
            thetaList = linspace(0,pi/2,nFrames+1); % rotation
            thetaList = thetaList(1:nFrames);
    end

    rawImage.cdata = uint8(zeros([imageMax 3]));
    rawImage.colormap = [];
    tic
    for iFrame = 1:nFrames
        
        
        disp([num2str(iFrame) ' of ' num2str(nFrames)])
        image = zeros(imageMax);
        
        %% calculate perspective
        
        
        theta = thetaList(iFrame);
        R = [cos(theta) -sin(theta);
            sin(theta) cos(theta)];
        x_y_rot = R * x_y; % rotate
        
        plotImCoordX = (x_y_rot(1,:)./depth*imageMax(2))+imageMax(2)/2; % plane of view at depth 1;  focal point at [0 0 0]
        plotImCoordY = (x_y_rot(2,:)./depth*imageMax(1))+imageMax(1)/2; % plane of view at depth 1;  focal point at [0 0 0]
        
        distEuclid = sqrt( x_y_rot(1,:).^2 + x_y_rot(2,:).^2 + depth.^2 );%  euclidean distance focal point at [0 0 0]
        
        %% Rendering
        % It could be done with figure and plot functions but they are not
        % designed to display objects on the pixel and sub-pixel scale.
        % The result just looks bad.
        
        % So a own small rendering engine is used:
        % check if a particular star rendering is necessary
        % lay a grid over the star
        % find if grid-points are within the star circle
        % if yes, match the grid-point to a pixel and add its portion of
        % "brightness" to the image
        
        startIndex = find(depth>0,1,'first'); % depth was sorted
        stopIndex = find(depth<renderDist,1,'last'); % depth was sorted
        
        for k = startIndex:stopIndex % for all stars that would be in range
            
            radiusStar = sizze(k)/distEuclid(k)/2;
            
            if plotImCoordX(k)+radiusStar < 0  % outside to the left
                continue % next star
            end
            if plotImCoordX(k)-radiusStar > imageMax(2)+1  % outside to the right
                continue % next star
            end
            
            if plotImCoordY(k)+radiusStar < 0  % outside to the bottom
                continue % next star
            end
            if plotImCoordY(k)-radiusStar > imageMax(1)+1  % outside to the top
                continue % next star
            end
            
            if radiusStar*2>starRenderLimInPixel % diameter at least a defined part of a pixel
                window = ceil(radiusStar *subsampPix); % size of grid to map over stars
                
                inCircle = false(2*window,2*window);
                
                for ix = 1:window % grid-points in positive x
                    for iy = 1:window % grid-points in positive y
                        
                        distEuclidPixel = sqrt( ((ix-0.5)/subsampPix).^2 + ((iy-0.5)/subsampPix).^2);%  euclidian distance focal point at [0 0 0]
                        
                        % do a quarter of the circle and flip to -x,y  and -x,-y  and  x,-y (1/8 would be enough actually)
                        if  radiusStar>distEuclidPixel % is grid in reach of star
                            inCircle(window+ix  ,window+iy  ) = true; % 1 quadrant
                            inCircle(window+ix  ,window+1-iy) = true; % 2 quadrant
                            inCircle(window+1-ix,window+iy  ) = true; % 3 quadrant
                            inCircle(window+1-ix,window+1-iy) = true; % 4 quadrant
                            
                        end
                    end
                end
                xGrid = round((( (1:2*window)-0.5-window) /subsampPix)+plotImCoordX(k)); % create x-grid and round to nearest Image pixel
                yGrid = round((( (1:2*window)-0.5-window) /subsampPix)+plotImCoordY(k)); % create y-grid and round to nearest Image pixel
                
                for ix = 1:2*window % pixel grid in x
                    if and(xGrid(ix) >= 1, xGrid(ix) <= imageMax(2)) % if in image
                        for iy = 1:2*window % pixel grid in y
                            if and(yGrid(iy) >= 1, yGrid(iy) <= imageMax(1)) % if in image
                                if inCircle(ix,iy) % if in Circle
                                    image(yGrid(iy),xGrid(ix)) = image(yGrid(iy),xGrid(ix))+1/subsampPix^2; % add brightness to pixel
                                end
                            end
                        end
                    end
                end
            end
        end
        
        image(image>1) = 1; % fix "over exposure"
        
        depth = depth-speed; % move
        
        %% save animation
        f = rawImage;
        f.cdata(:,:,1) = uint8(image*255);
        f.cdata(:,:,2) = uint8(image*255);
        f.cdata(:,:,3) = uint8(image*255);
        
        map = gray(8); % self created map. Otherwise use; [im,map] = rgb2ind(f.cdata,4,'nodither')
        
        imtemp = rgb2ind(f.cdata,map,'nodither');
        im(:,:,1,iFrame) = imtemp;
        if iFrame == 1
            im(1,1,1,nFrames) = 0; % allocate
        end
        
    end
    
    toc
    imwrite(im,map,fullfile(pathstr, [saveName '.gif']),'DelayTime',1/30,'LoopCount',inf) %
    disp([saveName '.gif  has ' num2str(numel(im)/10^6 ,4) ' Megapixels']) % Category:Animated GIF files exceeding the 50 MP limit
end
rng('shuffle'); % reset rand