#!/usr/bin/perl

#
# File:    acf.pl
# Descr:   A perl script for creating an AUTO constants file.
# Author:  Warren Weckesser, Department of Mathematics, Colgate University
# Date:    15 June 2003
#
# Last update: 18 June 2003
#
print "AUTO Constants File Builder   (18 June 2003)\n\n";
print "This perl script creates an AUTO constants file.\n";
print "It assumes that you are interested in continuing equilibria,\n";
print "periodic orbits, or solutions of boundary value problems.\n";
print "\n";

###################################################

$prob = "";
while ($prob ne "a" && $prob ne "b" && $prob ne "c" && $prob ne "d")
   {
   print("Continuation problem to be solved:\n");
   print(" (a) Equilibria\n");
   print(" (b) Periodic solutions and rotations\n");
   print(" (c) Basic boundary value problem\n");
   print(" (d) Use the BVP solver to compute a Poincare map\n");
   $prob = &PromptInput(": ");
   }

###################################################
print("\nProblem constants...\n");

$NDIM = &PromptInput("NDIM: Dimension of the system: ");

if ($prob eq "a" || $prob eq "b")
   {
   $NBC = 0;
   $NINT = 0;
   }
elsif ($prob eq "c")
   {
   # Basic BVP. It can have a combination of boundary conditions
   # and integral conditions.
   #
   $NBCdef = $NDIM;
   $NBC = &PromptInput("NBC:  Number of boundary conditions: [$NBCdef] ",$NBCdef);
   $NINTdef = $NDIM - $NBC;
   $NINT = &PromptInput("NINT: Number of integral conditions: [$NINTdef] ",$NINTdef);
   if (($NBC + $NINT - $NDIM) != 0)
      {
      print("\nWARNING: NBC+NINT-NDIM is not zero.\nThis script is not designed to handle this case.  Continue at your own risk!\n");
      }   
   }
else # ($prob eq "d")
   {
   # Using the BVP solver to compute a Poincare map.
   # NDIM "boundary" conditions are actually the initial condition, and one
   # more boundary condition is the stopping condition.
   $NBC = $NDIM+1;
   $NINT = 0;
   }

###################################################
print("\nDiscretization constants...\n");

if ($prob eq "a")
   {
   # Don't need NTST and NCOL for equilibria (???)
   $NTST = 0;
   $NCOL = 0;
   }
else
   {
   # Periodic solutions or BVP, so we need discretization info.
   $NTSTdef = 25;
   $NTST = &PromptInput("NTST: Num. of mesh intervals used for discretization: [$NTSTdef] ",$NTSTdef);
   $NCOLdef = 4;
   $NCOL = &PromptInput("NCOL: Num. of Gauss colloc. points per mesh interval (2 <= NCOL <= 7): [$NCOLdef] ",$NCOLdef);
   }

###################################################
print("\nTolerances...\n");

$EPSLdef = 1.0e-6;
$EPSL = &PromptInput("EPSL: Newton/Chord rel. convergence criterion for equation parameters: [$EPSLdef] ",$EPSLdef);

$EPSUdef = 1.0e-6;
$EPSU = &PromptInput("EPSU: Newton/Chord rel. convergence criterion for solution components: [$EPSUdef] ",$EPSUdef);

$EPSSdef = 1.0e-4;
$EPSS = &PromptInput("EPSS: Rel. arclength convergence criterion for detecting special solutions: [$EPSSdef] ",$EPSSdef);

$ITMXdef = 8;
$ITMX = &PromptInput("ITMX: Max. number of iterations allowed in the location of special solutions (bifs, folds, etc): [$ITMXdef] ",$ITMXdef);

$ITNWdef = 5;
$ITNW = &PromptInput("ITNW: Max. number of Newton/Chord iterations: [$ITNWdef] ",$ITNWdef);

###################################################
print("\nStep size parameters...\n");

$DS = &PromptInput("DS:   Initial step size for the pseudo-arclength continuation: ");

$DSMIN = &PromptInput("DSMIN: Minimum pseudo-arclength step size: ");

$DSMAX = &PromptInput("DSMAX: Maximum pseudo-arclength step size: ");

###################################################
print("\nDiagram Limits...\n");

$NMX = &PromptInput("NMX: Max. number of steps along any branch: ");

$RL0 = &PromptInput("RL0: Lower bound on the principal continuation parameter: ");

$RL1 = &PromptInput("RL1: Upper bound on the principal continuation parameter: ");

$A0 = &PromptInput("A0: Lower bound on the principal solution measure: ");

$A1 = &PromptInput("A1: Upper bound on the principal solution measure: ");

###################################################
print("\nFree parameters...\n");

$ICP = &PromptInput("ICP: Enter the index of the continuation parameter: ");

###################################################
print("\nComputation constants...\n");

$ILPdef = 0;
$ILP = &PromptInput("ILP: Detection of folds (0=no, 1=yes): [$ILPdef] ",$ILPdef);

$ISPdef = 1;
$ISP = &PromptInput("ISP: Controls the detection of branch points, etc. (0, 1, 2, or 3): [$ISPdef] ",$ISPdef);

$ISWdef = 1;
$ISW = &PromptInput("ISW: Controls switching at branch points for diff. eqs. (1, -1 or 2): [$ISWdef] ",$ISWdef);

$IRSdef = 0;
$IRS = &PromptInput("IRS: Label of the solution where the continuation is to be restarted\n      (0 if the starting point is coded in the C file): [$IRSdef] ",$IRSdef);

if ($prob eq "a")
   {
   $IPS = 1;
   }
elsif ($prob eq "b")
   {
   $IPS = 2;
   }
else
   {
   $IPS = 4;
   }


###################################################
print("\nOutput control...\n");

$NPRdef = 0;
$NPR = &PromptInput("NPR: Save solutions every NPR steps (NPR=0 means don't ouput solutions): [$NPRdef] ",$NPRdef);

$IPLTdef = 0;
$IPLT = &PromptInput("IPLT: Chooses the principal output measure.\n   IPLT=0: plot the L2 norm\n   0 < IPLT <= NDIM: Plot the max of the IPLT'th solution component\n   (see the manual for more options)\n   : [$IPLTdef] ",$IPLTdef);

$ans = "";
while ($ans ne "y" && $ans ne "n")
   {
   $ans = &PromptInput("Are there any specific values of the control parameter for which you want a solution saved? (y or n) [n]: ","n");
   }
if ($ans eq "y")
   {
   $valstr = &PromptInput("Enter the parameters values, separated by spaces: ");
   @values = split(/\s+/,$valstr);
   }
 

############################################################
# The following are the defaults for parameters not set up
# by this script.

# According to the manual, the following are "strongly recommended",
# so we'll these defaults:
#  IAD=3, NWTN=3, IADS=1
$IAD  = 3;
$NWTN = 3;
$IADS = 1;
# According to the manual, IID=2 is the recommended value.
$IID = 2;

# This script doesn't handle NTHU
$NTHU = 0;

# This script does not set up files that use these constants:
$MXBF = 0;
$JAC  = 0;

############################################################


print("\nThat's it!  Enter a file name, or just hit enter to abort.\n");

$fn = &PromptInput("Constants file name: [cancel] ","");

if ($fn ne "")
   {
   open(FN,">$fn") || die "Sorry, can not open $fn for output.\n";
   printf FN "%d %d %d %d                NDIM,IPS,IRS,ILP\n", $NDIM,$IPS,$IRS,$ILP;
   if ($prob eq "a" || $prob eq "c")
      {
      $NICP = 1;
      printf FN "%d %s                 NICP,(ICP(I),I=0,NICP-1)\n",$NICP,$ICP;
      }
   else # ($prob eq "b" || $prob eq "d")
      {
      $NICP = 2;
      printf FN "%d %s 10                NICP,(ICP(I),I=0,NICP-1)\n",$NICP,$ICP;
      }
   printf FN "%d %d %d %d %d %d %d %d    NTST,NCOL,IAD,ISP,ISW,IPLT,NBC,NINT\n",$NTST,$NCOL,$IAD,$ISP,$ISW,$IPLT,$NBC,$NINT;
   printf FN "%d %g %g %g %g           NMX,RL0,RL1,A0,A1\n",$NMX,$RL0,$RL1,$A0,$A1;
   printf FN "%d %d %d %d %d %d %d     NPR,MXBF,IID,ITMX,ITNW,NWTN,JAC\n",$NPR,$MXBF,$IID,$ITMX,$ITNW,$NWTN,$JAC;
   printf FN "%g %g %g         EPSL,EPSU,EPSS\n",$EPSL,$EPSU,$EPSS;
   printf FN "%g %g %g %d            DS,DSMIN,DSMAX,IADS\n",$DS,$DSMIN,$DSMAX,$IADS;
   if ($prob eq "a" || $prob eq "c" || $prob eq "d")
      {
      printf FN "0  NTHL\n";
      }
   else
      {
      printf FN "1             NTHL\n";
      printf FN "10 0.0\n";
      }
   printf FN "0  NTHU\n";
   if ($ans eq "y")
      {
      $NUZR = @values;
      printf FN "%d          NUZR,((I,UZR(I)),I=0,NUZR-1)\n",$NUZR;
      for ($i = 0; $i < $NUZR; $i++)
         {
         printf FN "%d %g\n",ICP,@values[$i];
         }
      }
   else
      {
      printf FN "0          NUZR,((I,UZR(I)),I=0,NUZR-1)\n",$NUZR;
      }
   close(FN);
   }

############################################################

#
# This subroutine does not let the user enter an empty string
# (unless a default string was provided as a second argument).
#
sub PromptInput
   {
   local($prompt,$default) = @_;
   local($q) = "";
   $askagain = 1;
   while ($askagain)
      {
      print "$prompt";
      $q = <STDIN>;
      chop($q);
      if ($q eq "" && defined($default))
         {
         $q = $default;
         $askagain = 0;
         }
      elsif ($q ne "")
         {
         $askagain = 0;
         }
      }
   $q;
   }