#define DATUM            "2009_07_07"
/****************************************************
This is the file 'spule.c.2009_07_07', 
genuine name:    'spule.c'.         

**************************************************
   First version: spule.c.2009_06_23
   The standard site for this file is: $HOME/c-gcc/

   cc -lm spule.c -Wall -o ~/bin/spule

   Adjustable parameters in the lines with: */      /*P*/         /*            
*************************************************
still missing:  7.5 T coil adatai 

*************************************************
magnetic field in a solenoid: 1 A/mm = 4 PI Oersted      

magnetic field in the axis of a current loop:
Bz=0.5*I*mu_0*R^2*(R^2+a^2)^-3/2 

magnetic field in the centre of a current loop:
Bz=0.5*I*mu_0*R^-1

mu_0 = 4 PI 10^-7 Vs/Am
*************************************************/
 
#define __NO_MATH_INLINES    /* makes slower but, may be, more precise */ 

/* for monitoring all errors:           int ERROR = 0;   
   for stopping at the first error:     int ERROR = 1;   */
int ERROR = 0;                                                 /*P*/

int TEST = 0;             /* active if 1; non-active if 0 */   /*P*/ 


/*for calculation with "long double" please uncomment the following lines:*/  
#define DOUBLE long double 
#define SQRT sqrtl
#define ATAN2 atan2l
#define LOG logl
#define FABS fabsl
#define LOG logl
/**/

/*for calculation with "double" please uncomment the following lines: 
#define DOUBLE double 
#define SQRT sqrt
#define ATAN2 atan2
#define LOG log
#define FABS fabs
#define LOG log
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include <fenv.h>

#define MAX_LEN 256      /* max. number of input characters */
#define EFFZERO 1e-320    /* my definition of zero min: 1e-320 */
#define LIM_EXP 70.      /* limit for the function exp */
#define PI 3.1415926535897932384626433832795028841971693993751

/* function prototypes: */
DOUBLE WWW(DOUBLE r, DOUBLE z, DOUBLE cf);
void F(DOUBLE cc);
void FF(DOUBLE w, int i);
void read_file (void); 
void write_file (void); 
void calculate (DOUBLE, DOUBLE);
void grafika (void);
void test (char *text);
void create_sample_file (void);
void help (void);
void panel (void);

void message (char *text, char *text1);
void error (char *text);
void terminate (char *text);
void wait (char *text);


/* global declarations: */
char *filename;        /* name of the coil set-up file*/
char gnuplotfile[MAX_LEN];
int order=0;     /* max number of the Landen iterations: usually 5 */       
DOUBLE XX,ZZ,R1,R2,Z1,Z2, FX, FZ, W, INTX, INTZ, NNN[100], MMM[100];
int N, I;
DOUBLE A,BR,BZ,B,Z0Z,X0X,MU0;
int II,NN,K,ITER,KO;
DOUBLE ZZ1[50],ZZ2[50],RR1[50],RR2[50],JJ[50];
DOUBLE x1,x2,yy1,y2,z1,z2;
int pixel=0;      /* number of points in the plot: usually  1000 */            
DOUBLE MU0=4.E-7*PI;

int errorcount = 0;     /* for the function error */
int gnuplotversion;      /* 0 if >= 4.0; >0 if <=3.7 */           

/**********************    function: main    ***************/
int main(int argc, char *argv[]) {
  char *kontr;

  if (argc == 2)
    {filename=argv[1]; 
      if (strcmp(filename, "sample")==0) create_sample_file ();
      else if (strcmp(filename, "help")==0) help (); }   

  if (argc != 8) {message("\n You have made a mistake in the input,",
                          "please try again using the format:\n");
    message (" ", 
	     "spule <file> <x1> <y1> <z1> <x2> <y2> <z2>\n");
    message ("<file> =","name of the command file");
    message ("<x1>, <y1>, <z1> =",
	     "coordinates of the one end of the line in mm");
    message ("<x2>, <y2>, <z2> =",
	     "coordinates of the other end of the line in mm\n");

    fprintf(stderr,"For a prototype command file please type 'spule sample'.\n");
    fprintf(stderr,
	    "For more information about this program please type 'spule help'.\n");
    terminate ("  input error  ");}
  
  filename = argv[1];
  
  x1 = 1.e-3*strtod(argv[2], &kontr);
  if(kontr == argv[2]) terminate(" wrong second parameter\n");
  
  yy1 = 1.e-3*strtod(argv[3], &kontr);
  if(kontr == argv[3]) terminate(" wrong third parameter\n");

  z1 = 1.e-3*strtod(argv[4], &kontr);
  if(kontr == argv[4]) terminate(" wrong fourth parameter\n");
  
  x2 = 1.e-3*strtod(argv[5], &kontr);
  if(kontr == argv[5]) terminate(" wrong fifth parameter\n");

  y2 = 1.e-3*strtod(argv[6], &kontr);
  if(kontr == argv[6]) terminate(" wrong sixth parameter\n");
  
  z2 = 1.e-3*strtod(argv[7], &kontr);
  if(kontr == argv[7]) terminate(" wrong seventh parameter\n");

  panel ();
  read_file();
      
  MMM[1]=1/SQRT(2.);
  NNN[1]=.5;
  ITER=order;
  for (I=1; I<=ITER; I++){
    II = I+1;
    MMM[II] = (MMM[I]+NNN[I])/2.;
    NNN[II] = SQRT(MMM[I]*NNN[I]);
    if (FABS(MMM[II]-NNN[II])< EFFZERO) goto AI;}
 AI:write_file();
  grafika ();
  error ("terminal"); 
  return 0;}
/**********************    function: main    ***************/

/**********************    function: read_file    ***************/
void read_file () {
  /* the geometry and current of the coils will be read in:  */
  int k=0;
  float f1,f2,f3,f4,f5;
  char File [MAX_LEN];
  char line [MAX_LEN+1];
  FILE *file; 

  sprintf(File, "%s", filename);
   
  printf("coil#    RR1[mm]   RR2[mm]   ZZ1[mm]   ZZ2[mm]   JJ[A/mm^2] \n");    
  file=fopen(File,"r");
  if (file==NULL)  terminate(" The desired file does not exist.  ");
  else 
    while (fgets (line, MAX_LEN, file) != NULL) {
      if (line[0]!='#') {
	if (pixel == 0){     sscanf (line, "%f",&f1); pixel = f1;}
	else if (order == 0){sscanf (line, "%f",&f1); order = f1;}
	else
	  {sscanf (line, "%f%f%f%f%f",&f1,&f2,&f3,&f4,&f5); 
	k=k+1;
	RR1[k]=1.e-3*f1;
	RR2[k]=1.e-3*f2;
	ZZ1[k]=1.e-3*f3;
	ZZ2[k]=1.e-3*f4;
	JJ[k] =1e6*f5;
	printf ("%4d  %10.1f%10.1f%10.1f%10.1f     %8.0f\n",
		k,f1,f2,f3,f4,f5);} } }
  printf ("\n");
  fclose(file);
  NN=k;
  return;}
/**********************    function: read_file    ***************/


/**********************    function: write_file    ***************/
void write_file () {
  /* transforms the (x,y,z) coordinates to the (r,z) coordinates,
     calculates the field values by the function calculate, 
     transforms the (BR,BZ) values to (Bx, By, Bz) values,
     writes the file .gnu which will be displayed 
     in the function grafika by the gnuplot */
  
  FILE *file; 
  int i;
  DOUBLE x,y,z,r,Bx,By;

  sprintf(gnuplotfile, "%s.gnu", filename);

  /* against compiler warning message: */
  file = fopen(gnuplotfile, "w"); fclose( file);  

  file = fopen(gnuplotfile, "w"); 
  fprintf (file, 
  "# i  x[mm] y[mm] z[mm] Bx[Gauss] By[Gauss] Bz[Gauss] B[Gauss]\n");
  printf ( 
  "   i   x[mm]   y[mm]   z[mm]:    Bx[Gauss]    By[Gauss]    Bz[Gauss]:     B[Gauss]\n");

  for (i=0; i<=pixel; i++) { 
    x=x1+(x2-x1)*i/(double)pixel;
    y=yy1+(y2-yy1)*i/(double)pixel;
    z=z1+(z2-z1)*i/(double)pixel;
    r = SQRT(x*x + y*y);  
    calculate (r,z);    
    if (r>EFFZERO) {Bx = BR * x / r;  By = BR * y / r;}
    else           {Bx = 0;           By = 0;         } 
    fprintf (file, "%3i %10.0Lf  %10.0Lf  %10.0Lf %12.4Lf  %12.4Lf %12.4Lf  %12.4Lf \n", 
	     i, 1e3*x, 1e3*y, 1e3*z, 1e4*Bx, 1e4*By, 1e4*BZ, 1e4*B); 
    if ((i==0)||(i==pixel/2)||(i==pixel)||(i==pixel/4)||(i==(pixel/4)*3))
      printf(
      "%4i %7.1Lf %7.1Lf %7.1Lf: %12.2Lf %12.2Lf %12.2Lf: %12.2Lf \n", 
      i, x*1000., y*1000., z*1000., Bx*1.e4, By*1.e4, BZ*1.e4, B*1.e4);}
  fclose( file);}
/**********************    function: write_file    ***************/


/***********************  function: grafika **********************/
      /* GNUPLOT HOWTO:               t '<titel>' w l 
         (in newer gnuplot versions:) lw <linewidth> lt <colour> 
         (in older gnuplot versions:) lw <linewidth> <colour> 
                                      w e == with errorbars */
void grafika(void){
  FILE *pipe;

  /* the gnuplot version will be asked: */
  pipe = popen("gnuplot", "w");
  fprintf(pipe, "set terminal png \n");
  gnuplotversion =  pclose(pipe);
  /*    if (gnuplotversion == 0) message(
	"A new gnuplot version is present.","\n");
	else message
	("An old gnuplot version is present: no png files will be available.",
	"\n");
  */	/* ready: gnuplot version */
  
  
  pipe = popen("gnuplot -geometry 900x400+10+10", "w");
  fprintf (pipe, "set grid \n");
  fprintf (pipe, 
  "set title '(%.1Lf,%.1Lf,%.1Lf)             <%s>             (%.1Lf,%.1Lf,%.1Lf)' \n", 
	   1.e3*x1, 1.e3*yy1,1.e3*z1,filename,1.e3*x2,1.e3*y2,1.e3*z2);
  fflush(pipe);
  if (gnuplotversion == 0)
    fprintf (pipe,    
	     "plot [:][:] '%s' using ($1):($5) t 'Bx' w l lw 2 lt 1, " 
	     "'%s' using ($1):($6) t 'By' w l lw 2 lt 3, "
	     "'%s' using ($1):($7) t 'Bz' w l lw 2 lt 4, "
 	     "'%s' using ($1):($8) t 'B ' w l lw 2 lt 2 \n",
	     gnuplotfile, gnuplotfile, gnuplotfile, gnuplotfile);
  else
    fprintf (pipe,   
	     "plot [:][:] '%s' using ($1):($5) t 'Bx' w l lw 2    1, " 
	     "'%s' using ($1):($6) t 'By' w l lw 2    3, "
	     "'%s' using ($1):($7) t 'Bz' w l lw 2    4, "
 	     "'%s' using ($1):($8) t 'B ' w l lw 2    2 \n",
	     gnuplotfile, gnuplotfile, gnuplotfile, gnuplotfile);
  fflush(pipe); 
  fprintf(stderr,"\n");
  fprintf(stderr,"plot window help:  h \n");
  fprintf(stderr,"            mark zoom region:  right mouse knob \n");
  fprintf(stderr,"            activate zoom :    left mouse knob \n");
  fprintf(stderr,"            unzoom:            u \n");
  fprintf(stderr,"            y logscale toggle: l \n\n");
 
  wait ("quit program   ");
  fprintf (pipe, "q  \n");        
  fflush(pipe);
  fclose (pipe);}
/*************************  function: grafika ***************************/


/**********************    function: calculate    ***************/
void calculate (DOUBLE r, DOUBLE z) {
  int i,j;
  
  XX=r; ZZ=z; BR=0.; BZ=0.;
  for (i=1; i<=NN; i++) {    
    Z1=ZZ1[i]-ZZ; Z2=ZZ2[i]-ZZ;  
    R1=RR1[i];    R2=RR2[i];        
    F(MMM[1]);    INTX=FX;          INTZ=FZ; 
    F(NNN[1]);    INTX=(INTX+FX)/2; INTZ=(INTZ+FZ)/2;  
    F(NNN[2]);    INTX=INTX+FX;     INTZ=INTZ+FZ; 
    for (j=3; j<=II-1; j++) FF(NNN[j],j-1);
    BR= BR+(MU0*JJ[i]/MMM[II])*INTX/16.;  
    BZ= BZ-(MU0*JJ[i]/MMM[II])*INTZ/16.;}
  B = SQRT(BR*BR+BZ*BZ);               }
/**********************    function: calculate    ***************/


/**********************    function: WWW    ***************/
DOUBLE WWW(DOUBLE r, DOUBLE z, DOUBLE cf){  
  DOUBLE erg;
  erg = (r*r-2.*r*XX*cf+XX*XX+z*z); 
  if (erg<0.) erg=0;
  erg = SQRT(erg);
  return(erg); }
/**********************    function: WWW    ***************/


/**********************    function: F    ***************/
void F(DOUBLE CC){
  DOUBLE WU,SF,CF,W,W1,W2,WA11,WA12,WA21,WA22,WB11,WB12,WB21,WB22;
  DOUBLE WC11,WC12,WC21,WC22,WD11,WD12,WD21,WD22;
  DOUBLE nenner, zahler;
 
  CF=SQRT(4.*CC*CC-1. );
  if (FABS(CC*CC-.5)>1e-10) SF=SQRT( 2.-4.*CC*CC); else SF=0.;    
 
  WA11=WWW(R1,Z1, CF); 
  WB11=WWW(R1,Z1, SF);
  WC11=WWW(R1,Z1,-CF);
  WD11=WWW(R1,Z1,-SF);   
  WA12=WWW(R1,Z2, CF);
  WB12=WWW(R1,Z2, SF);
  WC12=WWW(R1,Z2,-CF);
  WD12=WWW(R1,Z2,-SF);
  WA21=WWW(R2,Z1, CF);
  WB21=WWW(R2,Z1, SF);   
  WC21=WWW(R2,Z1,-CF);
  WD21=WWW(R2,Z1,-SF);
  WA22=WWW(R2,Z2, CF);
  WB22=WWW(R2,Z2, SF);
  WC22=WWW(R2,Z2,-CF);
  WD22=WWW(R2,Z2,-SF);   
  
  /*W =(WA11+R1-XX*CF)*(WA22+R2-XX*CF)/((WA12+R1-XX*CF)*(WA21+R2-XX*CF));*/
  zahler = (WA11+R1-XX*CF)*(WA22+R2-XX*CF);
  nenner = (WA12+R1-XX*CF)*(WA21+R2-XX*CF);
  if (FABS(nenner)<EFFZERO) {W=1.; test("1 ");} else W=zahler/nenner;
  
  /*W1=(WB11+R1-XX*SF)*(WB22+R2-XX*SF)/((WB12+R1-XX*SF)*(WB21+R2-XX*SF));*/
  zahler = (WB11+R1-XX*SF)*(WB22+R2-XX*SF);
  nenner = (WB12+R1-XX*SF)*(WB21+R2-XX*SF);
  if (FABS(nenner)<EFFZERO) {W1=1.; test("2 ");} else W1=zahler/nenner;

  /*W2=(WC11+R1+XX*CF)*(WC22+R2+XX*CF)/((WC12+R1+XX*CF)*(WC21+R2+XX*CF));*/
  zahler = (WC11+R1+XX*CF)*(WC22+R2+XX*CF);
  nenner = (WC12+R1+XX*CF)*(WC21+R2+XX*CF);
  if (FABS(nenner)<EFFZERO) {W2=1.;test("3 ");} else W2=zahler/nenner;

  W =W*W2;  
  
  /*W2=(WD11+R1+XX*SF)*(WD22+R2+XX*SF)/((WD12+R1+XX*SF)*(WD21+R2+XX*SF));*/
  zahler = (WD11+R1+XX*SF)*(WD22+R2+XX*SF);
  nenner = (WD12+R1+XX*SF)*(WD21+R2+XX*SF);      
  if (FABS(nenner)<EFFZERO) {W2=1.; test("4 ");} else W2=zahler/nenner;  

  W1=W1*W2;    
     
  /*W1=SF*SF*LOG(W1)+CF*CF*LOG(W);*/
  if (W1<EFFZERO) {W1=0.; test("5 ");} else W1=SF*SF*LOG(W1);
  if (W<EFFZERO) {; test("6 ");} else W1=W1+CF*CF*LOG(W);

  W =W1*XX+(WA11+WA22-WA12-WA21)*CF+(WB11+WB22-WB12-WB21)*SF-
    (WC11+WC22-WC12-WC21)*CF-(WD11+WD22-WD12-WD21)*SF;
  FX=W*CC;
  
  W =(WA11+Z1)*(WA22+Z2)*(WA12-Z2)*(WA21-Z1);
  WU=(WA11-Z1)*(WA22-Z2)*(WA12+Z2)*(WA21+Z1);
  
  if (FABS(WU)<EFFZERO) {W=1.; test("7 ");} else W =W/WU;  
  if (W<EFFZERO) {W1=0.; test("8 ");} else W1=LOG(W)*CF;                      
  W =(WB11+Z1)*(WB22+Z2)*(WB12-Z2)*(WB21-Z1);   
  WU=(WB11-Z1)*(WB22-Z2)*(WB12+Z2)*(WB21+Z1);    
  
  if (FABS(WU)<EFFZERO) {W=1.; test("9 ");} else W =W/WU;       
  if (W<EFFZERO) {; test("10 ");} else W1=W1+LOG(W)*SF;            

  zahler = (WC11+Z1)*(WC22+Z2)*(WC12-Z2)*(WC21-Z1);
  nenner = (WC11-Z1)*(WC22-Z2)*(WC12+Z2)*(WC21+Z1);         
  if (FABS(nenner)<EFFZERO) {W=1.; test("11 ");} else W=zahler/nenner;
  
  if (FABS(W)<EFFZERO){; test("12 ");} else W1=W1-LOG(W)*CF;        
  
  zahler = (WD11+Z1)*(WD22+Z2)*(WD12-Z2)*(WD21-Z1);
  nenner = (WD11-Z1)*(WD22-Z2)*(WD12+Z2)*(WD21+Z1);    
  if (FABS(nenner)<EFFZERO) {W=1.; test("13 ");} else W=zahler/nenner;

  /*W1=W1-LOG(W)*SF;*/ 
  if (W<EFFZERO) {; test("14 ");} else W1=W1-LOG(W)*SF;  
  FZ=W1*XX/2.;  
  
  /*W =(R1-XX*CF+WA11)*(R1-XX*SF+WB11)*(R1+XX*CF+WC11)*(R1+XX*SF+WD11)/
    ((R2-XX*CF+WA21)*(R2-XX*SF+WB21)*(R2+XX*CF+WC21)*(R2+XX*SF+WD21));*/
  zahler =(R1-XX*CF+WA11)*(R1-XX*SF+WB11)*(R1+XX*CF+WC11)*(R1+XX*SF+WD11);
  nenner =(R2-XX*CF+WA21)*(R2-XX*SF+WB21)*(R2+XX*CF+WC21)*(R2+XX*SF+WD21); 
  if (FABS(nenner)<EFFZERO) {W=1.; test("15 ");} else W=zahler/nenner;
  
  /*W1=(R2-XX*CF+WA22)*(R2-XX*SF+WB22)*(R2+XX*CF+WC22)*(R2+XX*SF+WD22)/
    ((R1-XX*CF+WA12)*(R1-XX*SF+WB12)*(R1+XX*CF+WC12)*(R1+XX*SF+WD12));*/
  zahler =(R2-XX*CF+WA22)*(R2-XX*SF+WB22)*(R2+XX*CF+WC22)*(R2+XX*SF+WD22);
  nenner =(R1-XX*CF+WA12)*(R1-XX*SF+WB12)*(R1+XX*CF+WC12)*(R1+XX*SF+WD12);
  if (FABS(nenner)<EFFZERO) {W1=1.; test("16 ");} else W1=zahler/nenner;
       
  /*FZ=FZ-Z1*LOG(W)-Z2*LOG(W1);*/   
  if (W<EFFZERO) {; test("17 ");} else FZ=FZ-Z1*LOG(W);
  if (W1<EFFZERO) {; test("18 ");} else FZ=FZ-Z2*LOG(W1);

  if (FABS (XX)<EFFZERO)  {W=0; /*test("19 ");*/} else {
      W1= ATAN2 (Z1*(R1-XX*CF),(XX*SF*WA11))+
	ATAN2 (Z2*(R2-XX*CF),(XX*SF*WA22))-
	ATAN2 (Z2*(R1-XX*CF),(XX*SF*WA12))-
	ATAN2 (Z1*(R2-XX*CF),(XX*SF*WA21)); 
  W=W1*SF; 
    W1=ATAN2 (Z1*(R1-XX*SF),(XX*CF*WB11))+
      ATAN2 (Z2*(R2-XX*SF),(XX*CF*WB22))-
      ATAN2 (Z2*(R1-XX*SF),(XX*CF*WB12))-
      ATAN2 (Z1*(R2-XX*SF),(XX*CF*WB21));  
  W=W+W1*CF; 
    W1= ATAN2 (Z1*(R1+XX*CF),(XX*SF*WC11))+
      ATAN2 (Z2*(R2+XX*CF),(XX*SF*WC22))-
      ATAN2 (Z2*(R1+XX*CF),(XX*SF*WC12))-
      ATAN2 (Z1*(R2+XX*CF),(XX*SF*WC21)); 
  W=W+W1*SF; 
    W1=ATAN2 (Z1*(R1+XX*SF),(XX*CF*WD11))+
      ATAN2 (Z2*(R2+XX*SF),(XX*CF*WD22))-
      ATAN2 (Z2*(R1+XX*SF),(XX*CF*WD12))-
      ATAN2 (Z1*(R2+XX*SF),(XX*CF*WD21)); 
    W=W+W1*CF; }
  FZ=FZ+W*XX;   
  FZ=FZ*CC;                      }
/**********************    function: F    ***************/


/**********************    function: FF    ***************/
void FF(DOUBLE w, int i){
  DOUBLE  w1;

  /* printf("FF: w=%24.20Lf  i=%i \n", w,i); */
  w1=SQRT( w*w-NNN[i]*NNN[i] ); 
  if (i > 2){FF(w+w1,i-1); FF(w-w1,i-1);}
  else { 
    F(w+w1); INTX=INTX+FX; INTZ=INTZ+FZ; 
    F(w-w1); INTX=INTX+FX; INTZ=INTZ+FZ; 
    INTX=INTX/2; INTZ=INTZ/2.;       }      }
/**********************    function: FF    ***************/

/**************** function: create sample file ********************/
void create_sample_file (){
  FILE *file;   
  
  file  = fopen("spule_sample", "w");
  fprintf(file, 
	  "# \n"   
	  "# \n"
	  "# \n"
"#   This is a command file for the program spule (spule.c).  \n"   
	  "# \n"
	  "#   No empty lines are allowed.\n"
"#   All comment lines must contain a # character at the very first place.\n"
"#   Number of points to be calculated for the graphics (pixel):\n"
	  "         1000\n"
	  "#\n"
"#   Maximum number of Landen iterations (order): \n"
	  "            5\n"
"#   Coil definitions (circular coils; axes are along z, current density j:\n"
	  "# \n"
"#   radius1 [mm]   radius2 [mm]    z1 [mm]    z2 [mm]    j [Amp/mm^2]\n"
"    115            185             -112        -32       130\n"
"    115            185               32        112       130\n"
"     70            106             -109        -24       130\n"
"     70            106               24        109       130\n"
	  "#\n"
	  "#\n"
	  "#\n"
	  "#\n");    
  fclose (file);  
  message (
  "\nA command sample file 'spule_sample' has been written or rewritten.", "");
  message ("For a demonstration of the program, please type \n",
     "                 spule spule_sample -1000 0 0 1000 0 0 ");
  message ("or \n",
     "                 spule spule_sample 0 0 -1000 0 0 1000 ");
  exit(1);} 
/**************** function: create sample file ********************/

/*********************** function: help ***************************/
void help (){
  panel ();

  fprintf(stderr,"Literature: W. Bartky, Rev. Mod. Phys. 10, (1938) 264 \n");
  fprintf(stderr,"agm = arithmetic-geometrical mean \n");
  fprintf(stderr,"\n");
  fprintf(stderr,"If the obtained curves are troubled, please choose another\n" 
	  "path in the vicinity of the desired one. Troubled sites can \n"
	  "be the coil surfaces or bodies. \n");
  fprintf(stderr,"\n"
	  "For this program following other programs must be installed:\n"
	  "    gcc gnuplot.\n"
	  "For a new version, please compile the C code by typing:\n"
	  "    cc -lm spule.c -Wall -o ~/bin/spule\n"); 
  fprintf(stderr,"\n"
	  "Author:  Janos Major (major@mf.mpg.de)\n"
	  "\n"
"---------------------------------------------------------\n"
"The program 'spule' is free software; you can redistribute it and/or\n"
"modify it under the terms of the GNU General Public License as published \n"
"by the Free Software Foundation; either version 2 of the License, or (at \n"
	  "your option) any later version.\n"
	  "\n"
"This program is distributed in the hope that it will be useful, but \n"
"WITHOUT ANY WARRANTY; without even the implied warranty of\n"
"MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU \n"
	  "General Public License for more details.\n"
	  "\n"
"You can obtain a copy of the GNU General Public License\n"
"from   http://www.gnu.org/copyleft/gpl.html; if not, write to \n"
"the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, \n"
	  "Boston, MA  02110-1301, USA.\n"
  "---------------------------------------------------------\n");
  exit(1);} 
/*********************** function: help ***************************/


/*********************** function: panel ***************************/
void panel (){
  message ("\n\n\n",""); 
  message( "  ______________________________________________________\n",
           "|                                                      |");
  message (" |                        SPULE                         |\n",
           "|                                                      |");
  message (" |                                                      |\n",
           "|      This program calculates the magnetic-field      |"); 
  message (" |   distribution of cylindrical concentric system of   |\n",
           "|   air-coils as defined in the command file, by the   |");
  message (" |     method of Landen transformation (agm method).    |\n",
           "|  The field values (x, y, and z components, modulus)  |");
  message (" |  along a straight line are calculated and displayed. |\n",
           "|                                                      |");
  message (" |       Max-Planck-Institut für Metallforschung        |\n",
           "|     Heisenbergstr. 3, D-70569 Stuttgart, Germany     |");
  message (" |                                                      |\n",
	   "|            (program version: " DATUM ")             |");
  message (" |______________________________________________________|",
           "\n\n"); } 
/*********************** function: panel ***************************/


/**********************    function: terminate    ***************/
void terminate (char *text){                   /* Text bei Fehler */
message("terminate:", text);
exit (1);}
/**********************    function: terminate    ***************/


/**********************      function: message    ***************/
void message (char *text1, char *text2){        /*Text beim Suchen*/
fprintf(stderr, "%s %s \n", text1, text2);}
/**********************      function: message    ***************/


/**********************      function: error    ***************/
void error (char *ref){                /* reports on runtime errors 
                                          Usage:       error ("label");
                          #include <fenv.h> and #include <errno.h> and 
                          int errorcount = 0;  (global) are necessary.*/
#ifdef _FENV_H   
  int raised;  
#endif /* fenv.h*/
  char text[MAX_LEN];  

  strcpy(text, "errno error message (");
  strcat(text, ref);  strcat(text, ")");
  perror(text);

#ifdef _FENV_H 
  raised = fetestexcept 
    (FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID);
  /*(FE_INEXACT | FE_DIVBYZERO | FE_UNDERFLOW | FE_OVERFLOW | FE_INVALID);*/
  if (raised) {errorcount = errorcount + 1; 
  sprintf(text, "%-3i", errorcount);
  strcat(text, "fp error message (");
  strcat(text, ref);  strcat(text, "):");
  /*if (raised & FE_INEXACT)   strcat(text, " INEXACT ");*/  
  if (raised & FE_DIVBYZERO) strcat(text, " DIVBYZERO ");
  if (raised & FE_UNDERFLOW) strcat(text, " UNDERFLOW ");
  if (raised & FE_OVERFLOW)  strcat(text, " OVERFLOW ");
  if (raised & FE_INVALID)   strcat(text, " INVALID ");
  message (text, " ");    
  if (ERROR) terminate ("           Floating point error detected.");}
  /* if (raised) terminate("floating point error"); */
  feclearexcept (FE_ALL_EXCEPT);    /**/  
#endif /* fenv.h*/
  errno = 0; /**/                  }
/**********************      function: error    ***************/


/**********************     wait      ********************/
void wait (char *text){
  char ch[MAX_LEN];  
  fprintf(stderr, "Further with RETURN:  %s", text);
  fgets(ch, MAX_LEN-1, stdin);}
/**********************     wait      ********************/

/**********************     test      ********************/
void test (char *text){  
  if (TEST) fprintf(stderr, "%s ", text);}
/**********************     test      ********************/



/**********************   whistle      ********************/
void whistle (void){putchar('\a'); fflush(stdout);}
/**********************   whistle     ********************/