
FourierTrigSeries := module() option package, load=ModuleLoad, unload=ModuleUnload;  
local ModuleLoad, ModuleUnload, createEmptySeries, setSeriesCoefficients, setSeriesGeneralCoefficient, setSeriesIndex, setSeriesBounds, setSeriesOrthogonalSystem, setSeriesPeriod, getSeriesVariable, getSeriesIndex, getSeriesPeriod, getSeriesBounds, getSeriesGeneralCoefficient, getSeriesOrthogonalSystem, getSeriesCoefficients, getSeriesGenre, getSeriesOrthoPoly;
export FOURIERSERIES, SERIESORTHOGONALSYSTEM, Degree, Create, Coefficients, Copy, Add, Derivate, ConvertToSum, Truncate, ScalarMultiply, ChangeBasis, Evaluate, GetFourierSeries, SimplifyCoefficients, ExploreFourierSeriesCoefficients, GetPartialSum, DrawPartialSum, DrawPeriodicExtension;
global `print/FOURIERSERIES`, `type/FOURIERSERIES`, `type/SERIESORTHOGONALSYSTEM`;

ModuleLoad := proc()   # define the types for FOURIERSERIES and SERIESORTHOGONALSYSTEM
  `type/FOURIERSERIES`:=t->evalb(op(0,t)='FOURIERSERIES');
  `type/SERIESORTHOGONALSYSTEM`:=t->evalb(t::{'CosSinTrigP'(name,name), 'CosTrigP'(name,name), 'SinTrigP'(name,name), 'ExpTrigP2'(name,name)});
  `print/FOURIERSERIES`:=proc(S);      # pre-printing conversion
    ConvertToSum('FOURIERSERIES'(S)):
   end proc:
end proc;

ModuleUnload := proc()   # define the types for FOURIERSERIES and SERIESORTHOGONALSYSTEM
  `type/FOURIERSERIES`:='`type/FOURIERSERIES`';
  `type/SERIESORTHOGONALSYSTEM`:='`type/SERIESORTHOGONALSYSTEM`';
  `print/FOURIERSERIES`:='`print/FOURIERSERIES`';
end proc;


createEmptySeries:=proc() local S, i;    # creates empty series of type FOURIERSERIES
  S:=array(0 .. 1,[]);     # the structure is similar to OrthogonalSeries
  S[0]:=table(sparse,[("dim")=1,("table")=(table(sparse,[])),("general")=[0,0],("bounds")=[0, infinity]]);
  S[1]:=table([("variable")=NULL,("genre")='CosSinTrigP',("index")=NULL,("period")=2*Pi]);
  'FOURIERSERIES'(S);
end proc:

setSeriesCoefficients:=proc(S::FourierTrigSeries:-FOURIERSERIES, coeffs::list) local coeff, t;  
  t:=op(1, S)[0];
  t["table"]:=table(sparse, coeffs);
  NULL;
end proc:

setSeriesGeneralCoefficient:=proc(S::FourierTrigSeries:-FOURIERSERIES, coeff) local gen, t, c, genre;
  t:=op(1,S)[0];
  genre:=op(1,S)[1]["genre"];
  if not type(coeff, 'list') then
    if genre='CosTrigP' then
      c:=[coeff, 0]
    elif genre='SinTrigP' then
      c:=[0, coeff]
    else
      c:=[coeff, coeff]
    end if;
  else
    c:=coeff;
  end if;
  t["general"]:=c;
  NULL;
end proc:

setSeriesBounds:=proc(S::FourierTrigSeries:-FOURIERSERIES, bounds) local t;
  t:=op(1, S)[0];
  t["bounds"]:=bounds;
  NULL;
end proc:

# sets the orthogonal system, variable and index
setSeriesOrthogonalSystem:=proc(S::FourierTrigSeries:-FOURIERSERIES, sys::FourierTrigSeries:-SERIESORTHOGONALSYSTEM) local sysname, t;  
  sysname:=op(0,sys);
  t:=op(1,S)[(1)];
  t["variable"]:=op(2,sys);    # set the variable, genre (system) and index
  t["genre"]:=sysname;
  t["index"]:=op(1,sys);
  NULL;
end:

setSeriesPeriod:=proc(S::FourierTrigSeries:-FOURIERSERIES, period::realcons) local t;
  t:=op(1, S)[1];
  t["period"]:=period;
  NULL;
end proc:

getSeriesVariable:=proc(S::FourierTrigSeries:-FOURIERSERIES) local v;
  v:=op(1,S)[(1)]["variable"];
  v;
end proc:

getSeriesIndex:=proc(S::FourierTrigSeries:-FOURIERSERIES) local i;
  i:=op(1,S)[(1)]["index"];
  i;
end proc:

getSeriesPeriod:=proc(S::FourierTrigSeries:-FOURIERSERIES) local p;
  p:=op(1, S)[1]["period"];
  p;
end proc:

getSeriesBounds:=proc(S::FourierTrigSeries:-FOURIERSERIES) local i;
  i:=op(1,S)[(0)]["bounds"];
  i;
end proc:

getSeriesGeneralCoefficient:=proc(S::FourierTrigSeries:-FOURIERSERIES) local c;
  c:=op(1,S)[0]["general"];
  if not ('type=list' in {args}) then
    if op(1,S)[(1)]["genre"]='CosTrigP' then
      c:=c[1];
    elif op(1,S)[(1)]["genre"]='SinTrigP' then
      c:=c[2];
    end if;
  end if;
  c;
end proc:

getSeriesCoefficients:=proc(S::FourierTrigSeries:-FOURIERSERIES) local t, c;   # returns the list of coefficients
  t:=eval(op(1,S)[(0)]["table"]);
  c:=op(2, op(t));
  c;
end proc:

# returns the system used (with the variable and the index) in the form "genre(index, variable)"
getSeriesOrthogonalSystem:=proc(S::FourierTrigSeries:-FOURIERSERIES) local sysname, index, variable, t;  
  t:=op(1,S)[(1)];
  variable:=t["variable"];
  sysname:=t["genre"];
  index:=t["index"];
  apply(sysname, index, variable);
end:

getSeriesGenre:=proc(S::FourierTrigSeries:-FOURIERSERIES) local genre, t;  
  t:=op(1,S)[(1)];
  genre:=t["genre"];
  genre;
end:

# creates one term multiplied by index (indeces for two-parts terms)
getSeriesOrthoPoly:=proc(S::FourierTrigSeries:-FOURIERSERIES) local genre, coeff, variable, index, pol, p, c;
  coeff:=`if`(nargs=2, op(2,[args]), 1);
  genre:=op(1,S)[(1)]["genre"];
  variable:=op(1,S)[(1)]["variable"];
  index:=op(1,S)[(1)]["index"];
  p:=2*Pi/getSeriesPeriod(S);
  if not type(coeff, 'list') then
    if genre='CosTrigP' then
      c:=[coeff, 0]
    elif genre='SinTrigP' then
      c:=[0, coeff]
    else
      c:=[coeff, coeff]
    end if;
  else
    c:=coeff;
  end if;
  if genre in {'CosSinTrigP', 'SinTrigP', 'CosTrigP'} then
    pol:=c[1]*cos(p*index*variable)+c[2]*sin(p*index*variable);
  elif genre='ExpTrigP2' then
    pol:=c[1]*exp(I*p*index*variable)+c[2]*exp(-I*p*index*variable);
  else
    pol:=NULL;
  end if;
  pol;
end proc:



# returns the degree if trigonometric polynomial (could be non-integer). 
Degree:=proc(S::FourierTrigSeries:-FOURIERSERIES) local d, g, t, l, p;
  d:=0;
  p:=2*Pi/getSeriesPeriod(S);
  g:=getSeriesGeneralCoefficient(S,'type'='list');  # get the general term
  if g<>[0,0] then # 
    d:=getSeriesBounds(S)[2]*p;
  else
    l:=getSeriesCoefficients(S);      # get the list of coefficients
    d:=max(op(map(x->lhs(x),l)))*p;
  end if;
  d;
end proc:

# creates the FOURIERSERIES series
Create:=proc(cf, system::FourierTrigSeries:-SERIESORTHOGONALSYSTEM) local coeff, period, S, elem, elem2, i, v, lb, clength;
  period:=`if`((nargs=3) and (op(3, [args])::realcons), op(3, [args]), 2*Pi);
  lb:=0;
  coeff:=`if`(cf::set, cf, {cf});
  S:=createEmptySeries();
  setSeriesOrthogonalSystem(S, system);
  for elem in coeff do
    if type(elem, list(`=`)) then 
      if op(0, system)='CosTrigP' then
        elem2:=map(x->`if`(rhs(x)::list, lhs(x)=rhs(x), lhs(x)=[rhs(x),0]), elem);
      elif op(0, system)='SinTrigP' then
        elem2:=map(x->`if`(rhs(x)::list, lhs(x)=rhs(x), lhs(x)=[0,rhs(x)]), elem);
      else
        elem2:=elem;
      end if; 
      setSeriesCoefficients(S, elem2);
      if lb=0 then lb:=max(op(map(x->op(1,x), elem)))+1; end if;
    elif type(elem, list) then
      if op(0, system)='SinTrigP' then
        elem2:=[seq(k=[0, elem[k]], k=1..nops(elem))];
        if lb=0 then lb:=nops(elem)+1;end if;
      elif op(0, system)='CosTrigP' then
        elem2:=[seq(k=[elem[k+1], 0], k=0..(nops(elem)-1))];
        if lb=0 then lb:=nops(elem);end if;
      else
        if lb=0 then lb:=nops(elem);end if;
        elem2:=[0=`if`(type(elem[1],'list'), elem[1], [elem[1], 0]), seq(k=`if`(type(elem[k+1],'list'), elem[k+1], [elem[k+1], elem[k+1]]), k=1..(nops(elem)-1))];
      end if;
      setSeriesCoefficients(S, elem2);
    elif type(elem, `=`) and (lhs(elem)='general') then
      setSeriesGeneralCoefficient(S, rhs(elem));
    elif type(elem, `=`) then
      setSeriesBounds(S, [op(1..2, rhs(elem))]);
      lb:=-1;
    else
      setSeriesGeneralCoefficient(S, elem);
    end;
  end do; 
  if (op(0, system)='SinTrigP') and (lb=0) then
    lb:=1;
  end if;
  if lb>0 then
  setSeriesBounds(S, [lb, infinity]);
  end if;
  setSeriesPeriod(S, period);
  v := getSeriesVariable(S); 
  if not type(v, name) then 
    error("wrong type of variable: must be name")
  end if;
  i := getSeriesIndex(S);
  if not type(i, name) then
    error("wrong index type: must be name")
  end if;
  S;  
end proc:

# returns the coefficient specified by index value
Coefficients:=proc(S::FourierTrigSeries:-FOURIERSERIES) local index, idx, c, subst, b;
  index:=`if`((nargs>1) and (op(2,[args])::integer), op(2, [args]), -1);
  if (nargs=1) or ((nargs=2) and ('type=list' in {args})) then
    return(getSeriesGeneralCoefficient(S));
  elif index>=0 then
    b:=getSeriesBounds(S);
    if index<b[1] then 
      c:=op(1,S)[(0)]["table"][(index)];
    elif index<=b[2] then
      subst:=op(1,S)[(1)]["index"]=index;
      c:=subs(subst, op(1,S)[(0)]["general"]);
    else
      c:=0;
    end if;
    if not (c::list) then 
      c:=[c,0];
    end if;
    if not('type=list' in {args}) then 
      if index=0 then
        c:=c[1];
      elif op(1,S)[(1)]["genre"]='CosTrigP' then
        c:=c[1];
      elif op(1,S)[(1)]["genre"]='SinTrigP' then
        c:=c[2];
      end if;
    end if;
    return(c);
  else
    error("index must be non-negative"); 
  end if;
end:

# converts FOURIERSERIES into ordinary series
ConvertToSum:=proc(S::FourierTrigSeries:-FOURIERSERIES) local sum1, c, i, b;
  sum1:=0;
  b:=getSeriesBounds(S);
  for i from 0 to (b[1]-1) do
    sum1:=sum1+subs(op(1,S)[(1)]["index"]=i, getSeriesOrthoPoly(S, Coefficients(S,i,'type'='list')));
  end do;
  sum1:=eval(sum1);
  if getSeriesGeneralCoefficient(S,'type'='list')<>[0,0] then
    sum1:=sum1+Sum(getSeriesOrthoPoly(S, Coefficients(S,'type'='list')), op(1,S)[(1)]["index"]=b[1]..b[2]);
  end if;
  sum1;
end proc:

Copy:=proc(S::FourierTrigSeries:-FOURIERSERIES) local S2;   # returns the copy of a series
  S2:=array(0 .. 1,[]);
  S2[0]:=copy(op(1,S)[0]);
  S2[1]:=copy(op(1,S)[1]);
  S2[0]["table"]:=copy(op(1,S)[0]["table"]);
  'FOURIERSERIES'(S2);
end proc:

Truncate:=proc(S::FourierTrigSeries:-FOURIERSERIES, k::integer) local S2, t, b, i;    # truncate series after k-th coefficient
  if k<0 then 
    error("index must be non-negative"); 
  end if;
  S2:=Copy(S);   # make a copy
  b:=getSeriesBounds(S);
  t:=op(1,S2)[0]["table"];
  if k<b[2] then
    if k>b[1] then   # just resize the bounds
      setSeriesBounds(S2, [b[1],k]);
    elif k=b[1] then    #convert the one-term sum into standalone term
      t[k]:=Coefficients(S2,k,'type'='list');
      setSeriesGeneralCoefficient(S2, [0,0]);
      setSeriesBounds(S2, [k+1, b[2]]);
    else     # delete redundant coefficients
      setSeriesGeneralCoefficient(S2, [0,0]);
      for i from k+1 to b[1] do
        t[i]:=0;
      end do;
    end if;
  end if;
  S2;
end proc:

# adds two series of the same type
Add:=proc(Series1::FourierTrigSeries:-FOURIERSERIES, Series2::FourierTrigSeries:-FOURIERSERIES) local i, S, S1, S2, b1, b2, lb, ub, g1, g2, g, d1, d2, l, t, c1, c2, p, genre1, genre2, mul1, mul2;
  if (getSeriesPeriod(Series1)<>getSeriesPeriod(Series2)) or    # test the period, the variable and the index name
     (getSeriesIndex(Series1)<>getSeriesIndex(Series2)) or
     (getSeriesVariable(Series1)<>getSeriesVariable(Series2)) then
    error("Series are not compatible");
  end if;
  mul1:=1;mul2:=1;
  if nargs>2 then
    mul1:=op(3, [args]);
    if has(mul1, getSeriesVariable(Series1)) then
      error(cat("wrong multiplier type: depends on ",getSeriesVariable(Series1)));
    elif has(a, getSeriesIndex(Series1)) then
      error(cat("wrong multiplier type: depends on ",getSeriesIndex(Series1)));
    end if;
  end if;
  if nargs>3 then
    mul2:=op(4, [args]);
    if has(mul2, getSeriesVariable(Series2)) then
      error(cat("wrong multiplier type: depends on ",getSeriesVariable(Series2)));
    elif has(a, getSeriesIndex(Series2)) then
      error(cat("wrong multiplier type: depends on ",getSeriesIndex(Series2)));
    end if;
  end if;
  if mul1<>1 then
    S1:=ScalarMultiply(Series1, mul1);
  else
    S1:=Series1;
  end if;
  if mul2<>1 then
    S2:=ScalarMultiply(Series2, mul2);
  else
    S2:=Series2;
  end if;
  genre1:=getSeriesGenre(S1);
  genre2:=getSeriesGenre(S2);
  if (genre1<>genre2) and ('ExpTrigP2' in {genre1, genre2}) then    # test genre compatibility
    error("Series are not compatible");
  end if;
  S:=createEmptySeries();   # new series for the sum
  if ({'CosTrigP', 'SinTrigP'} subset {genre1, genre2}) or     # genre transformation
     ('CosSinTrigP' in {genre1, genre2}) then
    setSeriesOrthogonalSystem(S, 'CosSinTrigP'(getSeriesIndex(S1),getSeriesVariable(S1)));  # set the orthogonal system
  else
    setSeriesOrthogonalSystem(S, getSeriesOrthogonalSystem(S1));  # set the orthogonal system
  end if;
  setSeriesPeriod(S, getSeriesPeriod(S1));     # copy period
  b1:=getSeriesBounds(S1);
  b2:=getSeriesBounds(S2);
  g1:=getSeriesGeneralCoefficient(S1,'type'='list');  # get general terms
  g2:=getSeriesGeneralCoefficient(S2,'type'='list');
  g:=[0,0];
  p:=2*Pi/getSeriesPeriod(S1);
  d1:=Degree(S1)/p;  # get degrees and translate them to integer numbers
  d2:=Degree(S2)/p; 
  if (g1=[0,0]) and (g2=[0,0]) then    # both general coefficients are zero
    lb:=max(d1, d2)+1;    # set lower and upper bound
    ub:=0; 
  elif (g1<>[0,0]) and (g2<>[0,0]) then    # both general coefficients are nonzero
    if b1[2]>b2[2] then
      ub:=b1[2];
      if b1[1]>b2[2] then
        lb:=b1[1]
      else
        lb:=b2[2]+1;
      end if;
    elif b1[2]=b2[2] then
      ub:=b1[2];
      if b1[1]>=b2[1] then
        lb:=b1[1];
      else
        lb:=b2[1];
      end if;
    else 
      ub:=b2[2];
      if b1[2]>=b2[1] then
        lb:=b1[2]+1
      else
        lb:=b2[1];
      end if;
    end if:
  elif g1<>[0,0] then   # first general coefficient is nonzero
    ub:=b1[2];
    if d2>=b1[1] then
      lb:=d2+1;
    else
      lb:=b1[1];
    end if;
  else              # second general coefficient is nonzero
    ub:=b2[2];
    if d1>=b2[1] then
      lb:=d1+1;
    else
      lb:=b2[1];
    end if;
  end if;
  if ub=lb then     # if the sum has only one term, the sum is not necessary
    ub:=0; 
    lb:=lb+1;
  end if;
  # now we have new bounds
  t:=op(1,S)[(0)]["table"]; 
  for i from 0 to (lb-1) do   # set the coefficients
    c1:=Coefficients(S1,i,'type'='list');
    c2:=Coefficients(S2,i,'type'='list');
    t[i]:=c1+c2;   # add coefficients
  end do;
  if ub>lb then    # set the bounds and general coefficient
    setSeriesBounds(S, [lb, ub]);
    g1:=`if`(d1<lb, [0,0], g1);   
    g2:=`if`(d2<lb, [0,0], g2);
    setSeriesGeneralCoefficient(S, eval(g1+g2));   # add general coefficients
  else
    setSeriesBounds(S, [lb, 'infinity']);          # correct series bounds
  end if;
  S;
end proc:

# return series multiplied by expression independent on series variable and index
ScalarMultiply:=proc(S::FourierTrigSeries:-FOURIERSERIES, a::algebraic) local i, S2, l, g;
  if has(a, getSeriesVariable(S)) then
    error(cat("wrong multiplier type: depends on ",getSeriesVariable(S)));
  elif has(a, getSeriesIndex(S)) then
    error(cat("wrong multiplier type: depends on ",getSeriesIndex(S)));
  end if;
  S2:=Copy(S);
  g:=getSeriesGeneralCoefficient(S2,'type'='list');
  if a=0 then
    setSeriesBounds(S2, [0,'infinity']);
    setSeriesGeneralCoefficient(S2, [0,0]);
    setSeriesCoefficients(S2, []);
  else
    setSeriesGeneralCoefficient(S2, [a*g[1],a*g[2]]);
    l:=getSeriesCoefficients(S2,'type'='list');
    l:=map(x->lhs(x)=[a*rhs(x)[1],a*rhs(x)[2]], l);
    setSeriesCoefficients(S2, l);
  end if;
  S2;
end proc:

Derivate:=proc(S::FourierTrigSeries:-FOURIERSERIES, diff_var::name) local S2, genre, coeff, c, i, elem, e, c_val, g, gen, index, variable, bounds, p;
  S2:=NULL;
  genre:=op(1,S)[(1)]["genre"];    # get the system
  index:=getSeriesIndex(S);
  variable:=getSeriesVariable(S);
  bounds:=getSeriesBounds(S);
  coeff:=NULL;
  p:=2*Pi/getSeriesPeriod(S);
  c:=getSeriesCoefficients(S,'type'='list');   
  for e in c do
    c_val:=rhs(e);
    i:=lhs(e);   # get the index value
    if diff_var=variable then   # differentiate the series with respect to it's variable
      if i<>0 then
        # set new values of the coefficient
        elem:=`if`(genre='ExpTrigP2', [I*p*i*c_val[1], -I*p*i*c_val[2]], [p*i*c_val[2], -p*i*c_val[1]]);  
        coeff:=coeff, (i=elem);
      end if;
    else  # differentiate the series with respect to different variable
      coeff:=coeff, (i=map(diff, c_val, diff_var));
    end if;
  end do;     # now we have all new coefficients
  g:=getSeriesGeneralCoefficient(S,'type'='list');
  if diff_var=variable then   # differentiate the series with respect to it's variable
    if genre='ExpTrigP2' then
      gen:=[I*p*index*g[1], -I*p*index*g[2]];   # set new values 
      # create new series
      S2:=Create({[coeff], 'general'=gen, k=bounds[1]..bounds[2]}, 'ExpTrigP2'(index, variable), getSeriesPeriod(S));
    elif genre='CosTrigP' then
      gen:=[0, -p*index*g[1]];   # set new values 
      S2:=Create({[coeff], 'general'=gen, k=bounds[1]..bounds[2]}, 'SinTrigP'(index, variable), getSeriesPeriod(S));  
    elif genre='SinTrigP' then
      gen:=[p*index*g[2], 0];   # set new values 
      S2:=Create({[coeff], 'general'=gen, k=bounds[1]..bounds[2]}, 'CosTrigP'(index, variable), getSeriesPeriod(S));  
    else
      gen:=[p*index*g[2], -p*index*g[1]];   # set new values 
      S2:=Create({[coeff], 'general'=gen, k=bounds[1]..bounds[2]}, 'CosSinTrigP'(index, variable), getSeriesPeriod(S));  
    end if;
  else # differentiate the series with respect to different variable
    gen:=map(diff, g, diff_var);   # set new values 
    S2:=Create({[coeff], 'general'=gen, k=bounds[1]..bounds[2]}, apply(genre, index, variable), getSeriesPeriod(S));  
  end if;
  S2;
end proc:

# change the orthogonal system
ChangeBasis:=proc(S::FourierTrigSeries:-FOURIERSERIES, sys::FourierTrigSeries:-SERIESORTHOGONALSYSTEM) local S2, variable, index, old_genre, new_genre, old_c, old_gc, new_c, new_gc, old_coeffs, new_coeffs, elem, bounds, period, sintrigseries, costrigseries;
  S2:=NULL;
  index:=op(1, sys);
  variable:=op(2, sys);
  bounds:=getSeriesBounds(S);
  period:=getSeriesPeriod(S);
  if (getSeriesIndex(S)<>index) or
     (getSeriesVariable(S)<>variable) then
    error("Series are not compatible");
  end if;
  old_genre:=getSeriesGenre(S);  
  new_genre:=op(0, sys);
  new_coeffs:=NULL;
  new_c:=[0,0];
  old_gc:=getSeriesGeneralCoefficient(S, 'type'='list');
  new_gc:=[0,0];
  if (old_genre in {'CosSinTrigP', 'CosTrigP', 'SinTrigP'}) and (new_genre='ExpTrigP2') then
    old_coeffs:=getSeriesCoefficients(S);
    for elem in old_coeffs do    # compute new coefficients
      old_c:=rhs(elem);
      if lhs(elem)=0 then
        new_c:=old_c;
      else
        new_c[1]:=1/2*(old_c[1]-I*old_c[2]);
        new_c[2]:=1/2*(old_c[1]+I*old_c[2]);
      end if;
      new_coeffs:=new_coeffs, lhs(elem)=new_c;
    end do;
    new_gc[1]:=1/2*(old_gc[1]-I*old_gc[2]);
    new_gc[2]:=1/2*(old_gc[1]+I*old_gc[2]);
    S2:=Create({[new_coeffs], 'general'=new_gc, k=bounds[1]..bounds[2]}, 'ExpTrigP2'(index, variable), period);
  elif (old_genre='ExpTrigP2') and (new_genre in {'CosSinTrigP', 'CosTrigP', 'SinTrigP'}) then
    sintrigseries:=true;
    costrigseries:=true;
    old_coeffs:=getSeriesCoefficients(S);
    for elem in old_coeffs do    # compute new coefficients
      old_c:=rhs(elem);
      if lhs(elem)=0 then
        new_c:=[old_c[1], 0];
      else
        new_c[1]:=old_c[1]+old_c[2];
        new_c[2]:=I*(old_c[1]-old_c[2]);
      end if;
      if (new_c[1]<>0) then sintrigseries:=false;end if;   # trying to recognize system SinTrigP or CosTrigP
      if (new_c[2]<>0) then costrigseries:=false;end if;
      new_coeffs:=new_coeffs, lhs(elem)=new_c;
    end do;
    new_gc[1]:=old_gc[1]+old_gc[2];
    new_gc[2]:=I*(old_gc[1]-old_gc[2]);
    if (new_gc[1]<>0) then sintrigseries:=false;end if;   # trying to recognize system SinTrigP or CosTrigP
    if (new_gc[2]<>0) then costrigseries:=false;end if;
    S2:=Create({[new_coeffs], 'general'=new_gc, k=bounds[1]..bounds[2]}, 'CosSinTrigP'(index, variable), period);
    if sintrigseries=true then
      setSeriesOrthogonalSystem(S2, 'SinTrigP'(index, variable));
    elif costrigseries=true then
      setSeriesOrthogonalSystem(S2, 'CosTrigP'(index, variable));
    end if;
    if (new_genre='SinTrigP') and (sintrigseries=false) then
      print(`Cannot convert to SinTrigP, converting to CosSinTrigP`)
    elif (new_genre='CosTrigP') and (costrigseries=false) then
      print(`Cannot convert to CosTrigP, converting to CosSinTrigP`)
    end if;
  else
    error("unsupported conversion");
  end if;
  S2;
end proc:

Evaluate:=proc(S::FourierTrigSeries:-FOURIERSERIES) local point, index, S2;
  point:=NULL;
  index:=NULL;
  # parameters processing
  if (nargs=2) and type(args[2], `=`) and (lhs(args[2])=getSeriesVariable(S)) then
    point:=rhs(args[2]);
  elif (nargs=2) and type(args[2], `=`) and (lhs(args[2])='trunc') and type(rhs(args[2]), 'nonnegint') then
    index:=rhs(args[2]);
  elif (nargs=3) and type(args[2], `=`) and (lhs(args[2])=getSeriesVariable(S)) and type(args[3], `=`) and (lhs(args[3])='trunc') and type(rhs(args[3]), 'integer') and  (rhs(args[3])>=0) then
    point:=rhs(args[2]);
    index:=rhs(args[3]);
  elif nargs>1 then
    error("wrong type of parametes");
  end if;
  if index<>NULL then   # truncate a series
    S2:=ConvertToSum(Truncate(S, index));
  else
    S2:=ConvertToSum(S);
  end if;
  if point<>NULL then
    S2:=subs(getSeriesVariable(S)=point, S2);
  end if;
  eval(value(S2));
end proc:


SimplifyCoefficients:=proc(S::FourierTrigSeries:-FOURIERSERIES, fcn::{identical(collect), identical(expand), identical(factor), identical(normal), identical(simplify)}) local S2, coef, gen; 
  if fcn = 'collect' and nargs < 3 then 
    error("Wrong parameters: with `collect` at least 3 parameters are expected"); 
  end if; 
  S2 := Copy(S);
  coef:=getSeriesCoefficients(S);
  coef:=map(fcn, coef, args[3..-1]);
  SetSeriesCoefficients(S2);
  gen:=getSeriesGeneralCoefficient(S, 'type'='list');
  gen:=map(fcn, gen, args[3..-1]);
  setSeriesCoefficients(S2, coef);
  setSeriesGeneralCoefficient(S2, gen);
  S2;
end proc;


#`+`:=overload(
#  [
#    proc(s1::FourierTrigSeries:-FOURIERSERIES, s2::FourierTrigSeries:-FOURIERSERIES) option overload;
#      Add(s1, s2);
#    end proc,
#    proc(s1::FourierTrigSeries:-FOURIERSERIES, c::realcons) option overload; local s2, l, i, coef;
#      s2:=Copy(s1);
#      l:=NULL;
#      for coef in getSeriesCoefficients(s1) do
#        if lhs(coef)<>0 then
#          l:=l, coef;
#        end if;
#      end do;
#      if Coefficients(s1,0, 'type'='list') = [0,0] then
#        l:=l, 0=[c,0];
#      else
#        l:=l, 0=[Coefficients(s1,0, 'type'='list')[1]+c,0];
#      end if;
#      setSeriesCoefficients(s2, [l]);
#      s2;
#    end proc,
#    proc(c::realcons, s1::FourierTrigSeries:-FOURIERSERIES) option overload;
#      `+`(s1,c);
#    end proc
#  ]
#):

#`-`:=overload(
#  [
#    proc(s1::FourierTrigSeries:-FOURIERSERIES) option overload;
#      ScalarMultiply(s1, -1);
#    end proc
#  ]
#):

#`*`:=overload(
#  [
#    proc(s::FourierTrigSeries:-FOURIERSERIES, c::realcons) option overload;
#      ScalarMultiply(s, c);
#    end proc,
#    proc(c::realcons, s::FourierTrigSeries:-FOURIERSERIES) option overload;
#      ScalarMultiply(s, c);
#    end proc
#  ]
#):

#`/`:=overload(
#  [
#    proc(s::FourierTrigSeries:-FOURIERSERIES, c::realcons) option overload;
#      ScalarMultiply(s, 1/c);
#    end proc
#  ]
#):

#Diff:=overload(
#  [
#    proc(s::FourierTrigSeries:-FOURIERSERIES, vars::name) option overload; local l, s2;
#      l:=[args][2..-1];
#      if not type(l, list(name)) then
#        error("wrong type of arguments");
#      end if;
#      if l[1]=getSeriesVariable(s) then
#        s2:=Derivate(s);
#        if nops(l)>1 then
#          return(diff(s2, op(2..-1,l)));
#        else
#          return(s2);
#        end if;
#      else
#        return(0);
#      end if;
#    end proc
#  ]
#):


# computes the Fourier Series on the given interval
GetFourierSeries:=proc(func::algebraic) local f, f2, variable, interval, a, b, a0, an, bn, L, lb, i, m, sol1, s, coeffs, exp_form, max, S, prop;
global n;
  # test if first argument is a function
  if (nargs in {2,3,4}) and type(args[2], 'name'='realcons'..'realcons') then
    variable:=op(1, args[2]);  # set the variable and function
    f:=func;
  else
    error("wrong type of arguments");
  end if;
  interval:=op(2, args[2]);   # set the interval
  if (evalf(op(1,interval))>evalf(op(2,interval))) then
    a:=op(2,interval);
    b:=op(1,interval);
  else
    a:=op(1,interval);
    b:=op(2,interval);
  end if;
  exp_form:=`if`('exp' in {args}, true, false);  # set the type of series
  if ({'odd', 'even'} in {args}) or ((('odd' in {args}) or ('even' in {args})) and (a<>0) and (b<>0)) then
    error("wrong type of arguments");
  end if;
  if ('odd' in {args}) then
    prop:='odd';
    f2:=-subs(variable=-variable, func);
    if a=0 then
      a:=-b;
      f:=piecewise(variable<0, f2, func);
    else
      b:=-a;
      f:=piecewise(variable<0, func, f2);
    end if;
  elif ('even' in {args}) then
    prop:='even';
    f2:=subs(variable=-variable, func);
    if a=0 then
      a:=-b;
      f:=piecewise(variable<0, f2, func);
    else
      b:=-a;
      f:=piecewise(variable<0, func, f2);
    end if;
  else
    prop:=NULL;
  end if;
  L:=abs(b-a):  # length of the interval
  # find the coefficients
  a0:=eval(1/L*int(f, variable=a..b)):
  assume('m'::realcons):  
  an:=2/L*int(f*cos(2*Pi/L*m*variable), variable=a..b);
  bn:=2/L*int(f*sin(2*Pi/L*m*variable), variable=a..b);
  # trying to find index values when coefficients are not defined
  sol1:={solve(denom(an)=0, m)} union {solve(denom(bn)=0, m)};
  max:=0;
  for s in sol1 do
    if (type(s,'integer') and (s>max)) then 
      max:=s;
    fi;
  od;
  n:='n';
  coeffs:=(0=[a0,0]);
  for i from 1 to max do     # compute first coefficients 
    an:=2/L*int(f*cos(2*Pi/L*i*variable), variable=a..b);
    bn:=2/L*int(f*sin(2*Pi/L*i*variable), variable=a..b);
    coeffs:=coeffs, i=[subs(m=n, an), subs(m=n, bn)];
  od;
  # compute again general coefficients
  assume(m::integer, m>max);
  an:=2/L*int(f*cos(2*Pi/L*m*variable), variable=a..b);
  bn:=2/L*int(f*sin(2*Pi/L*m*variable), variable=a..b);
  an:=subs(m=n, an);
  bn:=subs(m=n, bn);
  if prop='odd' then
    S:=Create({[coeffs], 'general'=[an,bn], k=(max+1)..'infinity'}, 'SinTrigP'('n',variable), L);
  elif prop='even' then
    S:=Create({[coeffs], 'general'=[an,bn], k=(max+1)..'infinity'}, 'CosTrigP'('n',variable), L);
  else
    S:=Create({[coeffs], 'general'=[an,bn], k=(max+1)..'infinity'}, 'CosSinTrigP'('n',variable), L);
  end if;
  if exp_form then
    S:=ChangeBasis(S, 'ExpTrigP2'('n',variable));
  end if;
  S;
end:


# procedure computes Fourier series coefficients for the given function
ExploreFourierSeriesCoefficients:=proc(func::algebraic) local f, f2, variable, interval, a, b, a0, an, bn, L, lb, i, m, sol1, s, coeffs, exp_form, max, S, prop, coeff_table, an_conds, bn_conds, cn_conds1, cn_conds2, steps, k;
global n;
  # test if first argument is a function
  if (nargs in {2,3,4,5}) and type(args[2], 'name'='realcons'..'realcons') then
    variable:=op(1, args[2]);  # set the variable and function
    f:=func;
  else
    error("wrong type of arguments");
  end if;
  steps:=1;
  if (nargs>=3) and type(args[3], 'posint') then
    steps:=args[3];
  end if;
  interval:=op(2, args[2]);   # set the interval
  if (evalf(op(1,interval))>evalf(op(2,interval))) then
    a:=op(2,range);
    b:=op(1,range);
  else
    a:=op(1,interval);
    b:=op(2,interval);
  end if;
  exp_form:=`if`('exp' in {args}, true, false);  # set the type of series
  if ({'odd', 'even'} in {args}) or ((('odd' in {args}) or ('even' in {args})) and (a<>0) and (b<>0)) then
    error("wrong type of arguments");
  end if;
  if ('odd' in {args}) then
    prop:='odd';
    f2:=-subs(variable=-variable, func);
    if a=0 then
      a:=-b;
      f:=piecewise(variable<0, f2, func);
    else
      b:=-a;
      f:=piecewise(variable<0, func, f2);
    end if;
  elif ('even' in {args}) then
    prop:='even';
    f2:=subs(variable=-variable, func);
    if a=0 then
      a:=-b;
      f:=piecewise(variable<0, f2, func);
    else
      b:=-a;
      f:=piecewise(variable<0, func, f2);
    end if;
  else
    prop:=NULL;
  end if;
  L:=abs(b-a):  # length of the interval
  # find the coefficients
  a0:=eval(1/L*int(f, variable=a..b)):
  assume('m'::realcons):  
  an:=2/L*int(f*cos(2*Pi/L*m*variable), variable=a..b);
  bn:=2/L*int(f*sin(2*Pi/L*m*variable), variable=a..b);
  # trying to find index values when coefficients are not defined
  sol1:={solve(denom(an)=0, m)} union {solve(denom(bn)=0, m)};
  max:=0;
  for s in sol1 do
    if (type(s,'integer') and (s>max)) then 
      max:=s;
    fi;
  od;
  n:='n';
  coeffs:=(0=[a0,0]);
  for i from 1 to max do     # compute first coefficients 
    an:=2/L*int(f*cos(2*Pi/L*i*variable), variable=a..b);
    bn:=2/L*int(f*sin(2*Pi/L*i*variable), variable=a..b);
    coeffs:=coeffs, i=[subs(m=n, an), subs(m=n, bn)];
  od;
  # compute again general coefficients
  assume(m::integer, m>max);
  an:=2/L*int(f*cos(2*Pi/L*m*variable), variable=a..b);
  bn:=2/L*int(f*sin(2*Pi/L*m*variable), variable=a..b);
  an:=simplify(subs(m=n, an));
  bn:=simplify(subs(m=n, bn));  
  coeff_table:=table([coeffs]);
  if not exp_form then  # sincos form of FS
    print('a[0]'/2+Sum('a[n]'*cos(2*Pi*n*variable/L)+'b[n]'*sin(2*Pi*n*variable/L),n=1..infinity));
    an_conds:=n=0,2*a0;
    bn_conds:=NULL;
    for i from 1 to max do
      an_conds:=an_conds, n=i, coeff_table[i][1];
      bn_conds:=bn_conds, n=i, coeff_table[i][2];
    end do;
    if steps>1 then
      for i from 0 to (steps-1) do
          an_conds:=an_conds, n=steps*k+i, subs(n=steps*k+i,an);
          bn_conds:=bn_conds, n=steps*k+i, subs(n=steps*k+i,bn);
      end do;
    else
      an_conds:=an_conds, an;
      bn_conds:=bn_conds, bn;
    end if;
    if nops([an_conds])>1 then
      an_conds:=simplify(piecewise(an_conds)) assuming k::integer;
      print('a[n]'=an_conds);
    else
      print('a[n]'=an);
    end if;
    if nops([bn_conds])>1 then
      bn_conds:=simplify(piecewise(bn_conds)) assuming k::integer;
      print('b[n]'=bn_conds);
    else
      print('b[n]'=bn);
    end if;
  else  # exp form of the FS
    print('c[0]'+Sum('c[n]'*exp(I*2*Pi*n*variable/L)+'c[-n]'*exp(-I*2*Pi*n*variable/L),n=1..infinity));
    print('c[0]'=a0);
    cn_conds1:=NULL;
    cn_conds2:=NULL;
    for i from 1 to max do
      cn_conds1:=cn_conds1, n=i, 1/2(coeff_table[i][1]-I*coeff_table[i][2]);
      cn_conds2:=cn_conds2, n=i, 1/2(coeff_table[i][1]+I*coeff_table[i][2]);
    end do;
    if steps>1 then
      for i from 0 to (steps-1) do
          cn_conds1:=cn_conds1, n=steps*k+i, 1/2*(subs(n=steps*k+i,an)-I*subs(n=steps*k+i,bn));
          cn_conds2:=cn_conds2, n=steps*k+i, 1/2*(subs(n=steps*k+i,an)+I*subs(n=steps*k+i,bn));
      end do;
    else
      cn_conds1:=cn_conds1, 1/2*(an-I*bn);
      cn_conds2:=cn_conds2, 1/2*(an+I*bn);
    end if;
    if nops([cn_conds1])>1 then
      cn_conds1:=simplify(piecewise(cn_conds1)) assuming k::integer;
      print('c[n]'=cn_conds1);
    else
      print('c[n]'=1/2*(an-I*bn));
    end if;
    if nops([cn_conds2])>1 then
      cn_conds2:=simplify(piecewise(cn_conds2)) assuming k::integer;
      print('c[-n]'=cn_conds2);
    else
      print('c[-n]'=1/2*(an+I*bn));
    end if;
  end if;
end:


GetPartialSum:=proc(S::FourierTrigSeries:-FOURIERSERIES, index::nonnegint);
  Evaluate(S, trunc=index);
end proc:


DrawPartialSum:=proc(S::FourierTrigSeries:-FOURIERSERIES, index::nonnegint) local f;
  f:=unapply(Evaluate(S, trunc=index), getSeriesVariable(S));
  plot(f, args[3..-1]);
end proc:


DrawPeriodicExtension:=proc(func::algebraic, var_range::name=realcons..realcons, newrange::realcons..realcons) local a, b, f,f2, variable, r, prop, L, modL, g1, g2, a2, b2, d, d2, disc, limit1, limit2, i, j, k, pom1, pom2;
  variable:=op(1, var_range);
  f:=func;
  r:=op(2,var_range);
  if (evalf(op(1,r))>evalf(op(2,r))) then
    a:=op(2,r);
    b:=op(1,r);
  else
    a:=op(1,r);
    b:=op(2,r);
  end if;
  if ({'odd', 'even'} in {args}) or ((('odd' in {args}) or ('even' in {args})) and (a<>0) and (b<>0)) then
    error("wrong type of arguments");
  end if;
  if ('odd' in {args}) then
    prop:='odd';
    f2:=-subs(variable=-variable, func);
    if a=0 then
      a:=-b;
      f:=piecewise(variable<0, f2, func);
    else
      b:=-a;
      f:=piecewise(variable<0, func, f2);
    end if;
  elif ('even' in {args}) then
    prop:='even';
    f2:=subs(variable=-variable, func);
    if a=0 then
      a:=-b;
      f:=piecewise(variable<0, f2, func);
    else
      b:=-a;
      f:=piecewise(variable<0, func, f2);
    end if;
  else
    prop:=NULL;
  end if;
  L:=abs(b-a):  # length of the interval
  f:=unapply(f,variable);
  # we will draw limits at discont. points
  if 'drawlimits' in {args} then
    if evalf(op(1,newrange))<evalf(op(2,newrange)) then
      a2:=op(1,newrange);
      b2:=op(2,newrange);
    else
      a2:=op(2,newrange);
      b2:=op(1,newrange);
    end if;
    disc:=fdiscont(f(x),x=a..b);
    d:=NULL;
    d2:=NULL;
    limit1:=limit(f(x),x=a,right);
    limit2:=limit(f(x),x=b,left);
    if (not(type(limit1,infinity) or type(limit2,infinity))) and (limit1<>limit2) then
      d:=d,[a,evalf((limit1+limit2)/2)];
      d:=d,[b,evalf((limit1+limit2)/2)];
    fi;
    for i from 1 to nops(disc) do
      limit1:=limit(f(x),x=op(1,disc[i]),left);
      limit2:=limit(f(x),x=op(2,disc[i]),right);
      if (not(type(limit1,infinity) or type(limit2,infinity))) and (abs(limit1-limit2)>0.01) then
        d:=d,[op(1,disc[i]),evalf((limit1+limit2)/2)];
      fi;
    od;
    d:=[d];
    if (a<>0) then pom1:=floor(a2/a)*signum(a); else pom1:=floor(a2/L); fi;
    if (b<>0) then pom2:=ceil(b2/b); else pom2:=ceil(b2/L)*signum(b2); fi;
    for i from pom1 to pom2 by 1 do
      for j from 1 to nops(d) do
        k:=d[j][1]+i*L;
        if ((evalf(k)>=evalf(a2)) and (evalf(k)<=evalf(b2))) then
          d2:=d2,[k,d[j][2]];
        fi;  
      od;
    od;
  end if;
  modL:=x->x-floor((x-a)/L)*L:
  f2:=x->f(modL(x)):
  g1:=plot(f2, newrange, discont=true, op({args[4..-1]} minus {drawlimits,odd,even}));
  g2:=NULL;
  if 'drawlimits' in {args} then 
    for i from 1 by 2 to 9 do
      g2:=g2,plot([d2], style=[point], symbol=CIRCLE, symbolsize=i, op({args[4..-1]} minus {drawlimits,odd,even}));
    od;
  fi;
  plots[display](g1,g2)
end proc:


end module:


