function y = FSFun(x)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% y = FSFun(x)
% ____________________________________________________________________________
% 18.385 (Nonlinear Dynamical Systems).         MIT, Fall 1999. R. R. Rosales.
%
% User provided function for input into fourierSC or heatSln.        Arbitrary
% periodic functions F = F(x) --- of period 2*pi --- are allowed,   provided a
% file FSFun.m (with a function script) is present in the directory from which
% the user is running MatLab. The MatLab statement y = FSFun(x)  should return
% the values of F(x) in the array y [i.e. y(n)=F(x(n))],  for any array x with
% x(n) in [0, 2*pi].     The function in FSFun will be used if the appropriate
% option picked when running fourierSC or heatSln (a default is provided).
%                                                  ^^^^^^^^^^^^^^^^^^^^^
% A simple example of a FSFun.m file is:
%
% function y = FSFun(x)
% y = exp(-10*(x-pi).^2);
%
% Note the use of the vector operand .^, needed to get a vector valued answer.
% 
% NOTE:   the function whose Fourier series is being computed is ASSUMED to be
%         piecewise continuous.     The discontinuities must be separated by a
% distance of at least 2*pi/50 with the current resolution in the scripts. The
% number of them cannot be too large either.    The evaluation near a point of
% discontinuity --- say, near x_d --- must be  done in such a way that
%
%  (A) F(x) defined by continuous expressions for both x > x_d and x < x_d;
%      with some value assigned at x = x_d.
%  (B) The value assigned at x_d satisfies:
%
%             F(x_d) = Fl  or  F(x_d) = Fr or F(x_d) = 0.5*(Fl + Fr),
%
%      where: Fl and Fr are the limits of F(x) from the left and right at x_d.
%      These three alternatives are useful in preserving nice symmetry (odd or
%      even) properties of the function F when evaluating it on a discrete set
%      of points x(j).
%
% For example, a "square wave" may be evaluated by the script
%
% function y = FSFun(x)
% y = ones(1, length(x));
% index = find((x > 0.5) & (x < 2.5));
% y(index) = zeros(1, length(index));
%
% 1) The first statement here creates an array of ones of the same size as x.
% 2) The second finds the sets of indexes j such that 0.5 < x(j) < 2.5.
% 3) The third statement sets the value of y to zero at the j's fund in (2).
%
% A special value could be assigned at the discontinuity points.  For example,
% add the statements
% 
% Idis = find(x == 0.5);
% y(Idis) = CONSTANT;
%
% where CONSTANT = 0 or 0.5 or 1 (here F+ = 1 and F- = 0).  Note that Idis may
% be the empty array if there is no j for which x(j) = 0.5.    Then the second
% statement here does not affect y.
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('\n NOTE: the default choice for FSFun is being used.')
fprintf('\n Do help FSFun to see how to use a different function.\n')
%
% y = (x.^2).*(x-2*pi).^2;   % EXAMPLE with discontinuous 2nd derivatives.
%
% ------ Example with asymmetric square wave with discontinuity close to edge
%        of domain.   Used to test part of code that creates the refined grid
%        xRef in fourierSC.m
y = zeros(size(x));
n1 = min(10, ceil(length(x)/10));
n2 = ceil(length(x)/4);
y(n1:(n1+n2)) = ones(1, n2+1);
%
% EOF