/* program for thin plate using second method - dual nodes and data at centre 
   of elements & activity flags + f.d.

to run, type:

flin [-l<lambda> -w<sor w> -n<no. of iterations> -f( full )
                                                    <input file> <output file>] 

or just flin and use the menu */

#include <math.h>
#include <graph.h>
#include <stdio.h>
#include <signal.h>
#include "az:plane.h"
#include "az:general.h"
#include "az_fe:fe.h"

FILE *fp, *fq, *flog = NULL;

unsigned short fontbuf[FONTSIZE];	/* buffer for largest expected font */
char *fontname = "fmacs:font.visual";
int h, w;

char sol_flag;         /* if true solution drawn at each iteration */
char norm_flag;        /* true if norm calculated at each iteration */
char j_flag;           /* true if j calculated at each iteration */

char menu_flag;        /* true if main menu active */
char move_flag;        /* true if text line cleared when cursor moved */
char plate_flag;       /* true if thin plate model used */

float mu;
float mu_4th;
float lambda;
float lam_sqr;
float alpha;
float beta;
float p;
/* constants for iteration */
float mu24mi, mu24pl, mu12mi, mu12pl, mu8mi, mu8pl, mu4mi, mu4pl;

float gen_norm_lim;    /* general norm limit - for terminating iterations */
float fin_norm_lim;    /* final norm limit - for terminating iterations p = 0 */
float cur_norm_lim;    /* current norm limit - for terminating iterations */
int max_no_it;         /* maximun no. of iterations - for terminating iterat. */
int tot_no_iter;       /* total number of iterations */
int no_of_steps;       /* number of steps in p */

float w_sor;           /* sor parameter */
int grey;              /* current grey level for text & information regions */

int xl_menu;              /* left x for menu region */
int xr_menu;              /* right x for menu region */
int y_menu;               /* top y for menu region */
int nop;                  /* no. of options for menu */

float i0_vec[ M +1 ], cts_vec[ M +1 ], phi_vec[ M + 1 ];
float store_phi_vec[ M + 1 ];
float v_vec[ M + 1 ], g_vec[ M + 1 ], h_vec[ M + 1 ], d4_vec[ M + 1 ];
char act_flag_vec[ M + 1 ];

extern stop();

extern void calc_p_alpha_consts(), calc_p_beta_consts();
extern void begin();
extern double g_alpha(), g1_alpha();
extern double g_beta(), g1_beta();

extern void initialise();
extern void pla_initialise();
extern void pla_iterate();
extern void init_act_flags();

void supervise();


main(argc,argv)
char **argv;
int argc;
{
   begin(argc,argv);
}

void supervise( start_flag ) /* supervises iteration */
char start_flag;
{
   float *q, *r;
   double phi_norm;
   static float store_gen_norm_lim, store_fin_norm_lim;
   static float old_phi_vec[ M + 1 ];           /* stores phi to calc norm */
   static int ip, iter;
   int i;
   char exit_flag = FALSE;  /* used for exiting function when MB pressed */
   char cont_flag = FALSE;  /* used if continuing iteration after MB pressed */


   /* check no. of steps reasonable */

   if ( no_of_steps <= 0 )
   {
      putchar( BEEP );
      write_text( "no. of steps out of range" );
      return;
   }

   if ( start_flag )
   {
      /* start of iteration */

      if ( sol_flag ) /* needed if starting another iteration */
         clear_vector( P_GREEN, old_phi_vec, X0, Y0 );

      /* set initial values for phi_vec */

      write_text( "calculating initial values" );

      if ( plate_flag )
         pla_initialise();
      else
         initialise();

      copy_vector( phi_vec, old_phi_vec ); /* to clear phi if re-starting */

      tot_no_iter = 0;
      ip = no_of_steps + 1;

      if ( sol_flag )
      {
         draw_vector( P_GREEN, phi_vec, X0, Y0 );
         draw_key( P_GREEN );
      }
   }
   else
   {
      /* to continue at same p & iter */
      cont_flag = TRUE;
      exit_flag = FALSE; /* in case continuing after iter. finished */
      if ( sol_flag )
      {
         clear_vector( P_GREEN, old_phi_vec, X0, Y0 );
         draw_vector( P_GREEN, phi_vec, X0, Y0 );
         draw_key( P_GREEN );
      }
   }

   write_text( "iterating : use middle MB to interrupt" );

   book_keep( START, &phi_norm );

   /* vary p between 1 and 0, and iterate for each p */

   do
   {
      if ( !cont_flag )
      {
         --ip;
         /* set value of p */

         if ( ip == no_of_steps )
            p = 1.0;
         else
            p = (double)ip / (double)no_of_steps;

         /* calculate consts for current p */

         if ( plate_flag )
         {
            calc_p_alpha_consts();
            calc_p_beta_consts();
         }
         else
         {
            calc_pmem_consts();
         }

         clear_info(); /* clear prev info during pause for init and 1st iter */
         put_p(); /* write p to screen and log */

         /* initialise activity flags */

         init_act_flags();

         iter = 0;

         /* store norm limits - in case changed when continuing */

         store_gen_norm_lim = gen_norm_lim;
         store_fin_norm_lim = fin_norm_lim;
      }
      else
      {
         /* continuing at same p */

         cont_flag = FALSE;
     
         /* check if norm limits changed & initialise activity flags if
                                                              neccesary */

         if ( ( ( gen_norm_lim < store_gen_norm_lim ) && ( p != 0 ) ) ||
              ( ( fin_norm_lim < store_fin_norm_lim ) && ( p == 0 ) ) )
            init_act_flags();
      }

      /* set value of cur_norm_lim */

      if ( ip == 0 )
         cur_norm_lim = fin_norm_lim;
      else
         cur_norm_lim = gen_norm_lim;

      if ( sol_flag )
         copy_vector( phi_vec, old_phi_vec ); /* to clear phi */

      /* iterate for each p */

      do
      {
         iter++;
         tot_no_iter++;


         /* to calc norm */
         copy_vector( phi_vec, store_phi_vec );

         if ( plate_flag )
            pla_iterate();
         else
            mem_iterate();

         book_keep( iter, &phi_norm );

         if ( (mouse->b & 0377) == 2 )
            exit_flag = TRUE;

      }
      while ( ( phi_norm > cur_norm_lim ) && 
              ( iter < max_no_it ) && !exit_flag );

      if ( sol_flag )
      {
         clear_vector( P_GREEN, old_phi_vec, X0, Y0 );
         draw_vector( P_GREEN, phi_vec, X0, Y0 );
      }
   }      
   while ( ( ip > 0 ) && !exit_flag );

   if ( exit_flag ) 
   {
      write_text( "iteration interrupted" );
      copy_vector( phi_vec, old_phi_vec ); /* to clear phi if restarting */
   }
   else
   {
      write_text( "iteration finished" );
      book_keep( FINISH, &phi_norm );
      if ( !sol_flag ) /* so solution won't have been drawn */
      {
         draw_vector( P_GREEN, phi_vec, X0, Y0 );
         draw_key( P_GREEN );
      }

      copy_vector( phi_vec, old_phi_vec ); /* to clear phi if restarting */
   }
}