%% slopefieldf(S, fn1, n)
%% %% where S=dy/dx
%% %% fn1 = filename for gnmuplot
%% %% fn2 = filename for image
%% %% func = an additional function to be plotted (if you want one)
%% %% n = line segment size 


%%
%% Examples:
%% slopefield( "x^2", "myplot.gp", "myplot.png", "(0.333333)*x*x*x")
%% %% This produces the slopefield for y=(1/3)x^3+C
%% %% The slope segments length's are defaulted to 1 unit
%%
%% slopefield( "x^2", "myplot.gp", "myplot.png", "(0.333333)*x*x*x", 0.75 )
%% %% This produces the slopefield for y=(1/3)x^3+C
%% %% The slope segments length's are 0.75 units in length
%%
%%
%% slopefield( "x*cos(x)*sin(x)", "test.gp", "test.png", "0.125*( sin(2*x)-2*x*cos(2*x))", 0.75 )
%% 
%% mike pellegrino - 2015
%% Gnuplot Slopefield Generator
function slopefield(S, fn1="default.gp", fn2="default.png", func="13", segmentsize=1)
  myfile=fopen(fn1, "w");
  
  segmentsize/=2;
  fprintf( myfile, "reset\n");
  fprintf( myfile, "set samples 1000\n");
  fprintf( myfile, "set label 1 \"x\" at 12.25, 0.1\n");
  fprintf( myfile, "set label 2 \"y\" at -0.2, 12.4\n");
  fprintf( myfile, "unset key\n");
  fprintf( myfile, "set xrange [-12:12]\n");
  fprintf( myfile, "set yrange [-12:12]\n");
  fprintf( myfile, "set xtics 1\n");
  fprintf( myfile, "set ytics 1\n");
  fprintf( myfile, "set grid\n");
  fprintf( myfile, "set title right\n");
  fprintf( myfile, "set term png size 640,640\n");
  fprintf( myfile, "set output \"%s\"\n", fn2);
  fprintf( myfile, "set arrow from 0,0 to 0,12\n");
  fprintf( myfile, "set arrow from 0,0 to 0,-12\n");
  fprintf( myfile, "set arrow from 0,0 to 12,0\n");
  fprintf( myfile, "set arrow from 0,0 to -12,0\n");

  for i=1:23
    x=i-12;
    for j=1:23
      y=j-12;
      try
	m=eval(S);
	if m == Inf
	  m=realmax;
	endif
      catch
	fprintf( myfile, "# [ERR0R]\n" );
	m=realmax;
      end_try_catch
      M(i,j)=m;
    endfor
  endfor

  for i=1:23
    x=i-12;
    for j=1:23
      
      y=j-12;
      m=M(i,j);
      x1=x-cos(atan(m))*segmentsize;
      x2=x+cos(atan(m))*segmentsize;
      y1=y-sin(atan(m))*segmentsize;
      y2=y+sin(atan(m))*segmentsize;

      fprintf( myfile, "set arrow from %f,%f to %f,%f lc 6 lw 2\n", x1,y1,x2,y2);

    endfor
  endfor

  fprintf( myfile, "f(x)=%s\n", func);
  fprintf( myfile, "plot f(x) lc 8 lw 3\n");
  fclose( myfile );
endfunction
