Logo Search packages:      
Sourcecode: varkon version File versions  Download package

sur905.c

/********************************************************************/
/*                                                                  */
/*  This file is part of the VARKON Geometry Library.               */
/*  URL:  http://www.varkon.com                                     */
/*                                                                  */
/*  This library is free software; you can redistribute it and/or   */
/*  modify it under the terms of the GNU Library General Public     */
/*  License as published by the Free Software Foundation; either    */
/*  version 2 of the License, or (at your option) any later         */
/*  version.                                                        */
/*                                                                  */
/*  This library is distributed in the hope that it will be         */
/*  useful, but WITHOUT ANY WARRANTY; without even the implied      */
/*  warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR         */
/*  PURPOSE.  See the GNU Library General Public License for more   */
/*  details.                                                        */
/*                                                                  */
/*  You should have received a copy of the GNU Library General      */
/*  Public License along with this library; if not, write to the    */
/*  Free Software Foundation, Inc., 675 Mass Ave, Cambridge,        */
/*  MA 02139, USA.                                                  */
/*                                                                  */
/*  (C)Microform AB 1984-1999, Gunnar Liden, gunnar@microform.se    */
/*                                                                  */
/********************************************************************/

#include "../../DB/include/DB.h"
#include "../include/GE.h"

/*!******************************************************************/
/*                                                                  */
/*  Function: varkon_cur_develop                   File: sur905.c   */
/*  =============================================================   */
/*                                                                  */
/*  Purpose                                                         */
/*  -------                                                         */
/*                                                                  */
/*  The function finds a developable ruling from a given point on   */
/*  curve to a point on the other curve.                            */
/*                                                                  */
/*  Input data are two curves and two parameter intervals for       */
/*  the search of the developable ruling. One interval shall        */
/*  have zero length, which identifies the given point on the       */
/*  given curve. The other interval must for the moment be          */
/*  within one segment of the curve.                                */
/*                                                                  */
/*  Input is also a pointer to a structure which defines the        */
/*  calculation criterion for the points (the objective             */
/*  function), the method for the numerical calculation,            */
/*  the number of restarts, the tolerances, etc.                    */
/*                                                                  */
/*  Note that the calling function handles the segments             */
/*  in the curve. Input parameters should normally be               */
/*  in the same segment. Also curve start and end                   */
/*  parameters should normally be modified in order                 */
/*  to respect the extension of the curve distance                  */
/*  ctol (the identical points criterion). Refer                    */
/*  to function varkon_cur_intplane (sur720).                       */
/*                                                                  */
/*  The methods available are Newton-Rhapson, ....                  */
/*                                                                  */
/*  Author                                                          */
/*                                                                  */
/*  Gunnar Liden                                                    */
/*                                                                  */
/*  Revisions                                                       */
/*                                                                  */
/*  1999-01-07   Originally written                                 */
/*  1999-12-01   Free source code modifications                     */
/*                                                                  */
/*                                                                  */
/******************************************************************!*/


/* ------------- Short description of function -----------------*/
/*                                                              */
/*sdescr varkon_cur_develop    SUR_RULED: Find developable ruler*/
/*                                                              */
/*------------------------------------------------------------- */

/*!-------------- Function calls (internal) ------------------------*/
/*                                                                  */
static short fdevrul();       /* Fctn f_cur and dfdt for DEVELRUL   */
static short addsol();        /* Add a non-existent solution.       */
static short fanalys();       /* Analysis of DEVELRUL solutions     */
#ifdef  DEBUG
static short initial();       /* Initialization of variables        */
#endif
/*                                                                  */
/*-----------------------------------------------------------------!*/

/*!-- Static (common) variables for the functions in this file -----*/
/*                                                                  */
static DBVector r_p;              /* Point   from point def. curve  */
static DBVector t_p;              /* Tangent from point def. curve  */
static DBint   ctype;             /* Case of computation            */
static DBfloat c_zero;            /* The zero (0.0) end criterion   */
static DBint   method;            /* Method which shall be used     */
static DBfloat t_s,t_e;           /* Param. interv. for calcul.     */
static DBfloat tglobal;           /* Current  parameter value       */
static DBfloat f_cur;             /* Current (objective) fctn value */
static DBfloat f_pre;             /* Previous function value        */
static DBfloat dfdt;              /* Function derivative w.r.t t    */
static DBint   numsol;            /* Number of solutions (=*p_npts) */
static EVALC   xyz;               /* Coordinates and derivatives    */
                                  /* for a point on a curve         */
static EVALC   xyz_p;             /* For analysis of DEVELRUL solut.*/
static DBfloat ctol;              /* Coordinate end criterion       */
static DBfloat ntol;              /* Angle      end criterion       */
static DBfloat comptol;           /* Computer tolerance (accuracy)  */
static DBfloat idpoint;           /* Identical points criterion     */
static DBint   niter;             /* Number of iterations           */
static DBfloat u_a[50];           /* All U parameter values         */
static DBfloat f_a[50];           /* All objective function values  */
static DBfloat d_a[50];           /* All objective derivat. values  */
/*-----------------------------------------------------------------!*/

/*!-------------- Function calls (external) ------------------------*/
/*                                                                  */
/* GE109                * Curve evaluation routine                  */
/* varkon_ini_evalc     * Initialize EVALC                          */
/* varkon_erpush        * Error message to terminal                 */
/*                                                                  */
/*-----------------------------------------------------------------!*/

/*!------------ Error messages and warnings ------------------------*/
/*                                                                  */
/* SU2943 = Called function xxxxxx failed in varkon_cur_develop     */
/* SU2983 = sur905 Illegal computation case=   for varkon_cur_deve  */
/* SU2973 = Internal function xxxxxx failed in varkon_cur_develop   */
/* SU2993 = Program error in varkon_cur_develop (sur905). Report !  */
/*                                                                  */
/*-----------------------------------------------------------------!*/


/*!****************** Function **************************************/
/*                                                                  */
   DBstatus  varkon_cur_develop  (

/*-------------- Argument declarations -----------------------------*/
/*                                                                  */
/* In:                                                              */
   DBCurve *p_cur1,      /* Curve 1                           (ptr) */
   DBSeg   *p_seg1,      /* Coefficients for curve p_cur1     (ptr) */
   DBfloat  t_in1[2],    /* Parameter curve segment limit p_cur1    */
                         /* (0):   Start U    (1): End   U          */
   DBCurve *p_cur2,      /* Curve 2                           (ptr) */
   DBSeg   *p_seg2,      /* Coefficients for curve p_cur2     (ptr) */
   DBfloat  t_in2[2],    /* Parameter curve segment limit p_cur2    */
                         /* (0):   Start U    (1): End   U          */
   IRUNON  *p_comp,      /* Pointer to computation data             */
/* Out:                                                             */
   DBint   *p_npts,      /* Number of output points           (ptr) */
   EVALC    xyz_a[SMAX], /* Array with all solutions                */
   DBint    d_code[SMAX])/* Developability analysis code            */
                         /* Eq.  1: Developable in both directions  */
                         /* Eq.  2: Developable in +U, cone in -U   */
                         /* Eq.  3: Developable in -U, cone in +U   */
                         /* Eq.  4: Developable but cones in +/-U   */
                         /* Eq.  0: No developable ruling d_code[0] */
                         /* Eq. -1: Failure to classify developable */
/*                                                                  */
/*-----------------------------------------------------------------!*/

{ /* Start of function */

/*!--------------- Internal variables ------------------------------*/
/*                                                                  */
   DBCurve *p_cur;       /* Curve for the numerical solution  (ptr) */
   DBSeg   *p_seg;       /* Coefficients for curve p_cur      (ptr) */
   DBCurve *p_curp;      /* Curve defining r_p and t_p        (ptr) */
   DBSeg   *p_segp;      /* Coefficients for curve p_curp     (ptr) */
   DBint    nstart;      /* The number of restarts                  */
   DBint    maxiter;     /* Maximum number of iterations            */
   DBfloat  deltat;      /* The parameter t step for the restarts   */
   DBint    i_r;         /* Loop index for the restarts             */
   DBint    j_n;         /* Loop index for the numerical solution   */
   DBfloat  h;           /* Step determined by Newton Rhapson       */

/*                                                                  */
/*-----------------------------------------------------------------!*/

   char     errbuf[80];  /* String for error message fctn erpush    */
   DBint    status;      /* Error code from a called function       */
#ifdef  DEBUG
   DBint    ix1;         /* Temporarely used loop index             */
#endif

/*--------------end-of-declarations---------------------------------*/

/*!                                                                 */
/* Algorithm                                                        */
/* =========                                                        */
/*                                                                 !*/

#ifdef DEBUG
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 Enter t_in1  %f %f  t_in2 %f %f\n", 
   t_in1[0],t_in1[1],  t_in2[0],t_in2[1]);
  }
#endif

/*!                                                                 */
/* 1. Check of input data and initializations                       */
/* __________________________________________                       */

/*  Initialize output coordinates and derivatives for DEBUG on.     */
/*  Call of initial, which also prints out the input data.          */
/*                                                                 !*/

/*  Printout of input data is in function initial                   */

#ifdef DEBUG
    status=initial(t_in1,t_in2,p_comp,p_npts,xyz_a,d_code);
    if (status<0) 
        {
        sprintf(errbuf,"initial%%sur905");
        return(varkon_erpush("SU2973",errbuf));
        }
#endif

/*!                                                                 */
/* Case of computation, method and number of restarts               */
   ctype   = p_comp->dev_un.ctype;
   method  = p_comp->dev_un.method;
   nstart  = p_comp->dev_un.nstart;
   maxiter = p_comp->dev_un.maxiter;
/*                                                                 !*/

/*!                                                                 */
/* Retrieve computation data for the input computation case         */
/* and definition of the end criterion for c_zero                   */
/*                                                                 !*/

   if      ( ctype  == DEVELRUL )
       {
       ctol    = p_comp->dev_un.ctol;
       ntol    = p_comp->dev_un.ntol;
       comptol = p_comp->dev_un.comptol;
       idpoint = p_comp->dev_un.idpoint;
/*    End criterion for the numerical solution                       */
       c_zero  = ntol;
       c_zero  = 10.0;    /*  Temporary   */
       }
   else
       {
       sprintf(errbuf,"(ctype )%%sur905");
       return(varkon_erpush("SU2993",errbuf));
       }


/*!                                                                 */
/* Determine which curve it is that shall be used in computation    */
/* Parameter interval t_s and t_e for the numerical solution        */
/*                                                                 !*/
   if      ( fabs(t_in1[1]-t_in1[0]) < comptol )
     { 
     p_cur   = p_cur2;
     p_seg   = p_seg2;
     p_curp  = p_cur1;
     p_segp  = p_seg1;
     if ( t_in2[0] < t_in2[1] )
       {
       t_s = t_in2[0];
       t_e = t_in2[1];
       }
     else
       {
       t_s = t_in2[1];
       t_e = t_in2[0];
       }
       xyz.t_global =  t_in1[0]; 
     }
   else if ( fabs(t_in2[1]-t_in2[0]) < comptol )
     { 
     p_cur   = p_cur1;
     p_seg   = p_seg1;
     p_curp  = p_cur2;
     p_segp  = p_seg2;
     if ( t_in1[0] < t_in1[1] )
       {
       t_s = t_in1[0];
       t_e = t_in1[1];
       }
     else
       {
       t_s = t_in1[1];
       t_e = t_in1[0];
       }
       xyz.t_global =  t_in2[0]; 
     }
   else
     {
     sprintf(errbuf,"No pt defined%%sur905");
     return(varkon_erpush("SU2993",errbuf));
     }


/*!                                                                 */
/* Point and tangent from curve with zero interval.                 */
/* Call of varkon_GE109 (GE109).                                  */
/*                                                                 !*/
   xyz.evltyp        = EVC_D2R; 

   status=GE109 ((DBAny *)p_curp, p_segp, &xyz );
   if (status<0) 
     {
     sprintf(errbuf,"GE109%%sur905");
     return(varkon_erpush("SU2943",errbuf));
    }

   xyz_p = xyz;
   r_p   = xyz_p.r;
   t_p   = xyz_p.drdt;


/*!                                                                 */
/* The parameter step for the restarts                              */
   deltat  = (t_e-t_s-comptol)/((DBfloat)nstart-1.0);
/* The number of output points                                      */
   numsol  = 0;
/*                                                                 !*/

#ifdef DEBUG
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 Master crv pt %6.4f Slave crv interv. t_s %6.4f t_e %6.4f\n",
   xyz_p.t_global, t_s,t_e);
  }
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 t_s %f t_e %f deltat %f nstart= %d ==> End start %f\n",
   t_s,t_e,deltat,(int)nstart,t_s+((DBfloat)nstart-1.0)*deltat);
  }
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 t_e - (nstart-1)*deltat %25.15f\n",
   t_e - t_s-((DBfloat)nstart-1.0)*deltat);
  }
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 r_p %f %f %f\n", r_p.x_gm, r_p.y_gm, r_p.z_gm );
  }
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 t_p %f %f %f\n", t_p.x_gm, t_p.y_gm, t_p.z_gm );
  }
#endif


/*!                                                                 */
/* 2. Numerical solution for a one parameter function               */
/* __________________________________________________               */

/* Start loop: Restarts i_r=1,2, .. nstart                          */
/*                                                                 !*/

   for (i_r = 1; i_r <= nstart ; i_r = i_r+1)
     {
/*!                                                                 */
/*   New parameter start value                                      */
/*   tglobal = t_s + deltat*(i_r-1)                                 */
/*   Start loop: Numerical i_n=1,2, ..                              */
/*                                                                 !*/

/*   Initialize number of iterations                                */
     niter = 0;

     tglobal = t_s + deltat*((DBfloat)i_r-1.0);

     for(j_n = 1; j_n <= maxiter ; j_n = j_n+1)
       {
/*!                                                                 */
/*     Calculate coordinates and derivatives.                       */
/*     Call of varkon_GE109   (GE109)                             */
/*                                                                 !*/
/*     (Error SU2943 if called function fails)                      */
/*     One additional iteration                                     */
       niter = niter + 1;

       xyz.t_global =  tglobal; 

       status=GE109 ((DBAny *)p_cur, p_seg, &xyz );
     if (status<0) 
       {
       sprintf(errbuf,"GE109 (loop)%%sur905");
       return(varkon_erpush("SU2943",errbuf));
    }

#ifdef DEBUG
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 r_x %f r_y %f r_z %f tglobal %9.6f\n",
      xyz.r.x_gm,   xyz.r.y_gm,   xyz.r.z_gm, tglobal);
  }
#endif

/*!                                                                 */
/*     Calculation of function f_cur and derivative dfdt.           */
/*     Call of fdevrul for ctype = DEVELSUR                         */
/*     (Error SU2973 if internal function fails)                    */
/*                                                                 !*/

       f_pre       = f_cur;  /* Tillfalligt !!! */

       if      ( ctype  == DEVELRUL )
         {
         status=fdevrul();
         if (status<0) 
            {
            sprintf(errbuf,"fdevrul%%sur905");
            return(varkon_erpush("SU2973",errbuf));
            }
         }
#ifdef DEBUG  
       else
         {
         sprintf(errbuf,"(ctype 2)%%sur905");
         return(varkon_erpush("SU2993",errbuf));
         }
#endif         






/*                                                                  */
/*     Calculated data to u_a, f_a and d_a.                         */
/*                                                                  */
       if ( niter >= 50)
          {
          sprintf(errbuf,"(niter>50)%%sur905");
          return(varkon_erpush("SU2993",errbuf));
          }

       u_a[niter-1]= tglobal;
       f_a[niter-1]= f_cur;
       d_a[niter-1]= dfdt;

/*!                                                                 */
/*     Is function value f_cur <= c_zero ?                          */
/*                                                                 !*/

#ifdef DEBUG
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 U %f f_cur %f c_zero %f dfdt %f\n"
   , tglobal, f_cur, c_zero, dfdt );
  }

if ( dbglev(SURPAC) == 2 && fabs(f_cur) < 10*c_zero )
{
for (ix1 = 1; ix1 <= niter; ix1 = ix1+1)
{
if ( fabs(f_a[ix1-1]) < 10000*c_zero )
{
  fprintf(dbgfil(SURPAC),
  "sur905 iteration %d U %25.20f f_cur %f dfdt %f\n"
   , (int)ix1  , u_a[ix1-1], f_a[ix1-1],d_a[ix1-1]);
}
}
}
#endif

      if( fabs(f_cur) <= c_zero ) 
          {
/*!                                                                 */
/*     Yes, f_cur= 0 for curve parameter value tglobal:             */
/*          Increase the number of solutions numsol and add         */
/*          the solution (R*3 and U,V) to output array xyz_a        */
/*          if the parameter tglobal is between t_s and t_e         */
/*          and if the solution not exists (call of addsol).        */
/*                                                                 !*/
/*   ?????? Note ?: tglobal= t_s and t_e-TOL4 will always be        */
/*                  since the minimum number of restarts is 2       */
/*  Att ta hand om !!!! oblem to have the solution as a u val */
/*           Accepting +/- TOL4 will give the wrong segment   */
/*           adress                                           */

          if ( tglobal >= t_s-TOL4)    /* Parameter is greater or   */
            {                          /* equal interval start t_s  */
            if ( tglobal <= t_e+TOL4)  /* Parameter value is less   */
              {                        /* than interval end t_e     */
              addsol(xyz_a);           /* Add non-existent solution */
              }                        /* End tglobal <= t_e        */
            }                          /* End tglobal >= t_s        */
          j_n = maxiter;               /* Refer to Note 2 above     */
          }                            /* End if f_cur <= c_zero    */

/*!                                                                 */
/*     Is derivative dfdt too small ?                               */
/*                                                                 !*/
          if(  fabs(dfdt) <= TOL3 )    /* Derivative <= TOL3 ?      */
            {                          /*                           */
/*!                                                                 */
/*     Yes, dfdt is less or equal TOL3:                             */
            j_n = maxiter;             /* ==> end of loop j_n       */
            h = TOL6;                  /* ??? Newton only ???       */
/*                                                                 !*/
            }                          /* End dfdt<= TOL3           */
          else                         /* Derivative >  TOL3        */
            {                          /*                           */
            h = -f_cur/dfdt;           /* Newton step ...           */
            }                          /* End if: dfdt <= TOL3      */

          /* Solution outside t_s < tglobal < t_e ? */
              if (tglobal < t_s + 0.01*TOL4 )
              {
              if ( h < 0 ) j_n = maxiter;
              }

          if ( tglobal > t_e + 0.01*TOL4 )
              {
              if ( h > 0 ) j_n = maxiter;
              }
          if ( fabs(h) > 2.0 ) j_n = maxiter;

#ifdef DEBUG
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 fabs(tglobal-t_s)= %f 0.01*TOL4 %f h %f j_n %d\n",
      fabs(tglobal-t_s),0.01*TOL4 , h , (int)j_n );
  }
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 fabs(tglobal-t_e)= %f 0.01*TOL4 %f h %f j_n %d\n",
      fabs(tglobal-t_e),0.01*TOL4 , h , (int)j_n );
  }
#endif

          /* Compute the next parameter value    */
          tglobal = tglobal + h;
          if ( tglobal < t_s - 0.01*TOL4 ) 
              tglobal = t_s + 0.001*TOL4/(DBfloat)niter;
          if ( tglobal > t_e + 0.01*TOL4 ) 
              tglobal = t_e - 0.001*TOL4/(DBfloat)niter;

#ifdef DEBUG
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 tglobal= %f   i_r= %d  j_n= %d\n",
      tglobal, (int)i_r, (int)j_n );
  }
#endif
        if ( j_n > maxiter )
            {
            sprintf(errbuf,"(maxiter)%% varkon_cur_develop (sur905)");
            return(erpush("SU2993",errbuf));
            }

        }   /* end loop j_n : Numerical solution    */
   }       /* end loop i_r : Restarts with deltat  */

/*!                                                                 */
/* 3. Exit                                                          */
/* ________                                                         */
/*                                                                  */
/* Number of solutions to output variable                           */
   *p_npts = numsol; 
/*                                                                 !*/

/*!                                                                 */
/* Analayse developable ruling. Call internal fctn fanalys.         */
/*                                                                 !*/

       if      ( ctype  == DEVELRUL )
         {
         status=fanalys(xyz_a,d_code);
         if (status<0) 
            {
            sprintf(errbuf,"fanalys%%sur905");
            return(varkon_erpush("SU2973",errbuf));
            }
         }
#ifdef DEBUG  
       else
         {
         sprintf(errbuf,"(ctype 3)%%sur905");
         return(varkon_erpush("SU2993",errbuf));
         }
#endif         




#ifdef DEBUG
if ( dbglev(SURPAC) == 1 )
{
if (*p_npts> 0 )
{
for (ix1 = 1; ix1 <= *p_npts; ix1 = ix1+1)
{
fprintf(dbgfil(SURPAC),
"sur905 Point %d U %9.6f X %9.2f Y %9.2f Z %9.2f Code %d\n",
      (int)ix1 , xyz_a[ix1-1].t_global  ,
    xyz_a[ix1-1].r.x_gm,xyz_a[ix1-1].r.y_gm,xyz_a[ix1-1].r.z_gm,
     (int)d_code[ix1-1]);
}
}
fflush(dbgfil(SURPAC)); 
}
if ( dbglev(SURPAC) == 1 && 0 == *p_npts )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 No solution! Master pt %6.4f Slave interv. %6.4f %6.4f\n",
   xyz_p.t_global, t_s,t_e);
  fflush(dbgfil(SURPAC)); 
  }
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 Exit *p_npts= %d\n", (int)*p_npts );
  fflush(dbgfil(SURPAC)); 
  }
#endif

    return(SUCCED);

} /* End of function                                                */

/*!****************************************************************!*/

#ifdef DEBUG
/*!********* Internal ** function ***********************************/
/*                                                                  */
/* Purpose                                                          */
/* --------                                                         */
/*                                                                  */
/* The function initializes the output variables and the            */
/* static (common) variables to F_UNDEF and I_UNDEF.                */

   static short initial(t_in1,t_in2,p_comp,p_npts,xyz_a,d_code)
/*                                                                  */
/******************************************************************!*/


/*!------------- Argument declarations -----------------------------*/
/*                                                                  */
/* In:                                                              */
   DBfloat  t_in1[2];    /* Parameter curve segment limits          */
                         /* (0):   Start U    (1): End   U          */
   DBfloat  t_in2[2];    /* Parameter curve segment limits          */
                         /* (0):   Start U    (1): End   U          */
   IRUNON  *p_comp;      /* Pointer to computation data             */
   int     *p_npts;      /* Number of output points (a pointer)     */
   EVALC    xyz_a[];     /* Array with all solutions                */
   DBint    d_code[];    /* Developability analysis code            */

/* Out:                                                             */
/*-----------------------------------------------------------------!*/

{ /* Start of function */


/*!--------------- Internal variables ------------------------------*/
/*                                                                  */
/*-----------------------------------------------------------------!*/

   DBint  i_sol;         /* Loop index solution                     */
   char   errbuf[80];    /* String for error message fctn erpush    */

/*--------------end-of-declarations---------------------------------*/

/*!                                                                 */
/* Algorithm                                                        */
/* =========                                                        */

/* 1. Printout of input variables                                   */
/*                                                                 !*/

if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905 Enter varkon_cur_develop*initial ******\n");
  fprintf(dbgfil(SURPAC),
  "sur905 Start t_in1(0) %f End t_in1(1) %f\n",
   t_in1[0],t_in1[1]);
  fprintf(dbgfil(SURPAC),
  "sur905 Start t_in2(0) %f End t_in2(1) %f\n",
   t_in2[0],t_in2[1]);
  fprintf(dbgfil(SURPAC),
  "sur905 Case of comp. %d  Method %d  No restarts %d Max iter %d\n",
      (int)p_comp->dev_un.ctype,
      (int)p_comp->dev_un.method, 
      (int)p_comp->dev_un.nstart, 
      (int)p_comp->dev_un.maxiter);
  if ( p_comp->dev_un.ctype == DEVELRUL)
      {
      fprintf(dbgfil(SURPAC),
      "sur905 Case DEVELRUL:  \n");
      }
  }

 fflush(dbgfil(SURPAC));

/*! The number of restarts must be between 2 and maxiter            */
/*  Error SU2993 if not OK.                                        !*/
    if (p_comp->dev_un.nstart > 1 || 
        p_comp->dev_un.nstart <= p_comp->dev_un.maxiter ) 
        {
        ;
        }
    else
        {
        sprintf(errbuf,"(nstart)%%sur905*initial");
        return(varkon_erpush("SU2993",errbuf));
        }

/*!                                                                 */
/* 2. Initialize output variables and static (internal) variables   */
/*                                                                 !*/
    *p_npts     = I_UNDEF;      

    r_p.x_gm    =   F_UNDEF;
    r_p.y_gm    =   F_UNDEF;
    r_p.z_gm    =   F_UNDEF;
    t_p.x_gm    =   F_UNDEF;
    t_p.y_gm    =   F_UNDEF;
    t_p.z_gm    =   F_UNDEF;

/* Initialize variabel EVALC and d_code                             */
/* Call of varkon_ini_evalc (sur776).                               */
/*                                                                 !*/

   for (i_sol = 1; i_sol <= SMAX;  ++i_sol)
     {
     varkon_ini_evalc (&xyz_a[i_sol-1]); 
     d_code[i_sol-1] = I_UNDEF;
     }

    return(SUCCED);

} /* End of function                                                */
/*!****************************************************************!*/
#endif

/*!************** Internal ** function ******************************/
/*                                                                  */
/* Purpose                                                          */
/* --------                                                         */
/*                                                                  */
/*  The function calcultes the function f_cur and dfdt              */
/*  for calculation case:  DEVELSUR                                 */
/*                                                                  */
               static  short   fdevrul()
/*                                                                  */
/*                                                                  */
/******************************************************************!*/


/*!------------- Argument declarations -----------------------------*/
/*                                                                  */
/* In: Global (static) variables:                                   */
/*     Points and derivatives xyz from varkon_cur_eval              */
/*     Point r_p and tangent t_p.                                   */
/* Out: Function f_cur and derivative dfdt                          */
/*                                                                  */
/*-----------------------------------------------------------------!*/

{ /* Start of function */

/*!--------------- Internal variables ------------------------------*/
/*                                                                  */
   DBfloat  dist_che;    /* Distance between curve points           */
 
/*-----------------------------------------------------------------!*/

   char   errbuf[80];    /* String for error message fctn erpush    */
/*!--------------- Theory ------------------------------------------*/
/*                                                                  */
/*  The function calculates a the objective function value for the  */
/*  creation of a ruled, developable surface                        */
/*                                                                  */
/*  Notations                                                       */
/*  r_p         = Input point (start of ruling from master curve)   */
/*  t_p         = Input tangent (from master curve for r_p)         */
/*  u           = Parameter for the slave curve                     */
/*  r_u         = Point on slave curve for a given u                */
/*  drdu        = First  derivative on slave curve                  */
/*  d2rdu2      = Second derivative on slave curve                  */
/*  f_obj       = Objective function                                */
/*  dfdu_obj    = Derivative of objective function w.r.t. u         */
/*                                                                  */
/*  The criterion for a developable ruling is that the              */
/*  surface normal direction is the same for all points             */
/*  along the ruling.                                               */
/*                                                                  */
/*  The ruling direction (vector) is (r_u - r_p) and the            */
/*  tangent drdu gives the surface normal direction in              */
/*  point r_u as the cross product dr_du X (r_u - r_p).             */
/*                                                                  */
/*  Equal directions for surface normals in r_p and r_u             */
/*  means that the tangent t_p in r_p must be perpendicular         */
/*  to the surface normal in r_u, i.e. must be in the               */
/*  tangent plane of r_u.                                           */
/*                                                                  */
/*  The scalar product of the direction of the surface normal       */
/*  (the cross product vector) and the tangent t_p is the           */
/*  objective function. The scalar product is zero when             */
/*  the vectors are perpendicular. Note that the objective          */
/*  function will be positive and negative around zero, i.e.        */
/*  it is not the minimum point we are trying to find.              */
/*                                                                  */
/*   Surface normal (cross vector product) in r_u:                  */
/*    _                                    _                        */
/*   !      e_x          e_y        e_z     !                       */
/*   !    drdu.x       drdu.y     drdu.z    ! =                     */
/*   ! r_u.x-r_p.x  r_u.y-r_p.y r_u.z-r_p.z !                       */
/*   !_                                    _!                       */
/*                                                                  */
/*   drdu.y*(r_u.z-r_p.z) - drdu.z*(r_u.y-r_p.y)*e_x                */
/*   drdu.z*(r_u.x-r_p.x) - drdu.x*(r_u.z-r_p.z)*e_y                */
/*   drdu.x*(r_u.y-r_p.y) - drdu.y*(r_u.x-r_p.x)*e_z                */
/*                                                                  */
/*  The objective function is the scalar product between            */
/*  the surface normal in r_u and the tangent in r_p:               */
/*  f_obj       =   r_u_surface_normal(u)*t_p =                     */
/*   (drdu.y*(r_u.z-r_p.z)-drdu.z*(r_u.y-r_p.y))*t_p.x +            */
/*   (drdu.z*(r_u.x-r_p.x)-drdu.x*(r_u.z-r_p.z))*t_p.y +            */
/*   (drdu.x*(r_u.y-r_p.y)-drdu.y*(r_u.x-r_p.x))*t_p.z ;            */
/*                                                                  */
/*  The derivative with respect to parameter u of the slave         */
/*  curve                                                           */
/*  dfdu_obj   =  d r_u_surface_normal/du * t_p =                   */
/*   (d2rdu2.y*(r_u.z-r_p.z)  - d2rdu2.z*(r_u.y-r_p.y)+             */
/*    drdu.y  *(drdu.z-0.0  ) - drdu.z  *(drdu.y-0.0 ))*t_p.x +     */
/*   (d2rdu2.z*(r_u.x-r_p.x)  - d2rdu2.x*(r_u.z-r_p.z)+             */
/*    drdu.z  *(drdu.x-0.0  ) - drdu.x  *(drdu.z-0.0 ))*t_p.y +     */
/*   (d2rdu2.x*(r_u.y-r_p.y)  - d2rdu2.y*(r_u.x-r_p.x)+             */
/*    drdu.x  *(drdu.y-0.0  ) - drdu.y  *(drdu.x-0.0 ))*t_p.z ;     */
/*                                                                  */
/*  dfdu_obj   =  d r_u_surface_normal/du * t_p =                   */
/*   (d2rdu2.y*(r_u.z-r_p.z)  - d2rdu2.z*(r_u.y-r_p.y)+             */
/*    drdu.y  *   drdu.z      - drdu.z  *   drdu.y     ) *t_p.x +   */
/*   (d2rdu2.z*(r_u.x-r_p.x)  - d2rdu2.x*(r_u.z-r_p.z)+             */
/*    drdu.z  *   drdu.x      - drdu.x  *   drdu.z     ) *t_p.y +   */
/*   (d2rdu2.x*(r_u.y-r_p.y)  - d2rdu2.y*(r_u.x-r_p.x)+             */
/*    drdu.x  *   drdu.y      - drdu.y  *   drdu.x     ) *t_p.z ;   */
/*                                                                  */
/*-----------------------------------------------------------------!*/

/*------------end-of-declarations-----------------------------------*/


/*!                                                                 */
/* Algorithm                                                        */
/* =========                                                        */
/*                                                                 !*/

/* Check distance between curves (curve points)                     */
/* Replace with better error message for the user  TODO             */
  dist_che = sqrt((xyz.r.x_gm-r_p.x_gm)*(xyz.r.x_gm-r_p.x_gm)+
                  (xyz.r.x_gm-r_p.x_gm)*(xyz.r.x_gm-r_p.x_gm)+
                  (xyz.r.x_gm-r_p.x_gm)*(xyz.r.x_gm-r_p.x_gm));
  if ( dist_che < 10.0*idpoint)
       {
       sprintf(errbuf,"(equal curves)%%sur905*fdevrul");
       return(varkon_erpush("SU2993",errbuf));
       }



  f_cur        = 
   (xyz.drdt.y_gm*(xyz.r.z_gm-r_p.z_gm)-
    xyz.drdt.z_gm*(xyz.r.y_gm-r_p.y_gm))*t_p.x_gm +
   (xyz.drdt.z_gm*(xyz.r.x_gm-r_p.x_gm)-
    xyz.drdt.x_gm*(xyz.r.z_gm-r_p.z_gm))*t_p.y_gm +
   (xyz.drdt.x_gm*(xyz.r.y_gm-r_p.y_gm)-
    xyz.drdt.y_gm*(xyz.r.x_gm-r_p.x_gm))*t_p.z_gm ;

  dfdt        = 
   (xyz.d2rdt2.y_gm*(xyz.r.z_gm-r_p.z_gm)  - 
    xyz.d2rdt2.z_gm*(xyz.r.y_gm-r_p.y_gm)+
    xyz.drdt.y_gm  *   xyz.drdt.z_gm      - 
    xyz.drdt.z_gm  *   xyz.drdt.y_gm     ) *t_p.x_gm +
   (xyz.d2rdt2.z_gm*(xyz.r.x_gm-r_p.x_gm)  - 
    xyz.d2rdt2.x_gm*(xyz.r.z_gm-r_p.z_gm)+
    xyz.drdt.z_gm  *   xyz.drdt.x_gm      - 
    xyz.drdt.x_gm  *   xyz.drdt.z_gm     ) *t_p.y_gm +
   (xyz.d2rdt2.x_gm*(xyz.r.y_gm-r_p.y_gm)  - 
    xyz.d2rdt2.y_gm*(xyz.r.x_gm-r_p.x_gm)+
    xyz.drdt.x_gm  *   xyz.drdt.y_gm      - 
    xyz.drdt.y_gm  *   xyz.drdt.x_gm     ) *t_p.z_gm ;


#ifdef DEBUG
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905*fdevrul tglobal= %f f_cur= %f dfdt= %f\n"
   , tglobal, f_cur, dfdt );
  }
#endif

       return(SUCCED);


} /* End of function */

/********************************************************************/








/*!************** Internal ** function ******************************/
/*                                                                  */
/* Purpose                                                          */
/* --------                                                         */
/*                                                                  */
/*  Add a solution if it is a new solution (not already existing).  */
/*                                                                  */
/*  Note that the identical points criterion (ctol) is used for     */
/*  all types of interrogation (c_zero is not used).                */
/*                                                                  */
               static  short   addsol(xyz_a)
/*                                                                  */
/*                                                                  */
/******************************************************************!*/


/*!------------- Argument declarations -----------------------------*/
/*                                                                  */
/* In:                                                              */
   EVALC   xyz_a[SMAX];  /* Array with all solutions                */

/*     Global (static) variables:                                   */
/*     Points and derivatives xyz from varkon_sur_eval              */
/* Out: Function                                                    */
/*                                                                  */
/*-----------------------------------------------------------------!*/

{ /* Start of function */

/*!--------------- Internal variables ------------------------------*/
/*                                                                  */
   DBint    i_sol;       /* Loop index for previous solutions       */
 
/*-----------------------------------------------------------------!*/

   DBfloat  dist;        /* Distance between solution points (R*3)  */
   char     errbuf[80];  /* String for error message fctn erpush    */

/*------------end-of-declarations-----------------------------------*/


/*!                                                                 */
/* Algorithm                                                        */
/* =========                                                        */
/*                                                                  */


#ifdef DEBUG
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905*addsol ctol %f numsol %d \n"
   ,   ctol , (int)numsol );
  }
#endif

/*!                                                                 */
/* 1. Check if solution already exists                              */
/* __________________________________                               */

/*  If numsol = 0 add solution (goto adds).                         */
/*                                                                 !*/

    if ( numsol == 0 ) goto  adds;

/*!                                                                 */
/*  Start loop previous solutions i_sol= 1,2, .. ,numsol.           */
/*                                                                 !*/

   for (i_sol = 1; i_sol <= numsol ; i_sol = i_sol+1)
     {

/*!                                                                 */
/*    Solution exists if the R*3 distance between xyz and xyz_a     */
/*    is less than ctol.  Goto sex if this is the case.             */
/*                                                                 !*/

     dist = SQRT (
       (xyz.r_x-xyz_a[i_sol-1].r_x)*(xyz.r_x-xyz_a[i_sol-1].r_x) +
       (xyz.r_y-xyz_a[i_sol-1].r_y)*(xyz.r_y-xyz_a[i_sol-1].r_y) +
       (xyz.r_z-xyz_a[i_sol-1].r_z)*(xyz.r_z-xyz_a[i_sol-1].r_z) );

#ifdef DEBUG
if ( dbglev(SURPAC) == 2 &&  dist > ctol )
  {
  fprintf(dbgfil(SURPAC),
  "sur905*addsol u %10.5f v %10.5f dist %f ctol %f \n",
   xyz_a[i_sol-1].u , xyz_a[i_sol-1].v , dist  ,   ctol );
  }
#endif

     if ( dist <= ctol ) goto sex;


     }  /* End loop i_sol= 1,2,..., numsol                          */



/*!                                                                 */
/* 2. Add solution                                                  */
/* _______________                                                  */


adds:  /* Label: Add solution                                       */

/*  Increase number of solutions.                                   */
    numsol = numsol + 1;
/*  Check that numsol is less or equal SMAX                         */
/*                                                                 !*/
    if ( numsol > SMAX )
       {
       sprintf(errbuf,"(SMAX)%% varkon_cur_develop (sur905)");
       return(varkon_erpush("SU2993",errbuf));
       }

/*!                                                                 */
/*  Points and coordinates to xyz_a.                                */
/*                                                                 !*/

    xyz_a[numsol-1]=xyz;     /* The R*3 solution point    */

sex:  /*! Label sex: Solution already exists                       */

    return(SUCCED);


} /* End of function */

/*!****************************************************************!*/


/*!************** Internal ** function ******************************/
/*                                                                  */
/* Purpose                                                          */
/* --------                                                         */
/*                                                                  */
/*  Analysis of DEVELRUL solutions                                  */
/*                                                                  */
               static  short  fanalys(xyz_a, d_code)
/*                                                                  */
/*                                                                  */
/******************************************************************!*/


/*!------------- Argument declarations -----------------------------*/
/*                                                                  */
/* In:                                                              */
   EVALC    xyz_a[SMAX]; /* Array with all solutions                */
/*     Global (static) variables:                                   */
/*     xyz_p                                                        */
/* Out:                                                             */
   DBint    d_code[SMAX];/* Developability analysis code            */
/*                                                                  */
/*-----------------------------------------------------------------!*/

{ /* Start of function */

/*!--------------- Internal variables ------------------------------*/
/*                                                                  */
   DBint    i_sol;       /* Loop index solution                     */
#ifdef  DEBUG
   DBfloat  f_cur_sol;   /* Function value f_cur for output point   */
   DBfloat  dfdt_sol;    /* Derivative of  f_cur for output point   */
   DBfloat  dfdt_switch; /* Derivative when curves are switched     */
#endif
   DBfloat  eps;         /* For approximation of r_p(t+eps) point   */
   DBfloat  f_cur_p;     /* Function value f_cur for + eps          */
   DBfloat  dfdt_p;      /* Derivative of  f_cur for + eps          */
   DBfloat  f_cur_n;     /* Function value f_cur for - eps          */
   DBfloat  dfdt_n;      /* Derivative of  f_cur for - eps          */
 
/*-----------------------------------------------------------------!*/

   short  status;        /* Error code from a called function       */
   char   errbuf[80];    /* String for error message fctn erpush    */

/*------------end-of-declarations-----------------------------------*/


/*!                                                                 */
/* Algorithm                                                        */
/* =========                                                        */
/*                                                                  */


#ifdef DEBUG
if ( dbglev(SURPAC) == 2 )
  {
  fprintf(dbgfil(SURPAC),
  "sur905*fanalys numsol %d \n"
   ,   (int)numsol );
  }
#endif

/*!                                                                 */
/* 1. Return if there are no solutions                              */
/* ___________________________________                              */
/*                                                                 !*/

    if ( 0 == numsol ) 
      {
      d_code[0] = 0;
      return(SUCCED);
      }

/* Initialization of internal variables                             */

#ifdef  DEBUG
   f_cur_sol    =  F_UNDEF;
   dfdt_sol     =  F_UNDEF;
   dfdt_switch  =  F_UNDEF;
   eps          =  F_UNDEF;
   f_cur_p      =  F_UNDEF;
   dfdt_p       =  F_UNDEF;
   f_cur_n      =  F_UNDEF;
   dfdt_n       =  F_UNDEF;
#endif

/*!                                                                 */
/* 2. Analysis w.r.t. "interchanged" curves                         */
/* ________________________________________                         */
/*                                                                 !*/

/*!                                                                 */
/*  Start loop solutions i_sol= 1,2, .. ,numsol.                    */
/*                                                                 !*/

   for (i_sol = 1; i_sol <= numsol ; i_sol = i_sol+1)
     {

/*!                                                                 */
/*    Recalculate function and derivative for + eps.                */
/*    Call of internal function fdevrul                             */
/*                                                                 !*/

      eps = 0.00001;
      r_p = xyz_p.r;
      t_p = xyz_p.drdt;
      r_p.x_gm=xyz_p.r.x_gm+
           eps*xyz_p.drdt.x_gm+eps*eps/2*xyz_p.d2rdt2.x_gm;
      r_p.y_gm=xyz_p.r.y_gm+eps*xyz_p.drdt.y_gm+
           eps*eps/2*xyz_p.d2rdt2.y_gm;
      r_p.z_gm=xyz_p.r.z_gm+eps*xyz_p.drdt.z_gm+
           eps*eps/2*xyz_p.d2rdt2.z_gm;
      t_p.x_gm=xyz_p.drdt.x_gm+eps/2*xyz_p.d2rdt2.x_gm;
      t_p.y_gm=xyz_p.drdt.y_gm+eps/2*xyz_p.d2rdt2.y_gm;
      t_p.z_gm=xyz_p.drdt.z_gm+eps/2*xyz_p.d2rdt2.z_gm;
      xyz = xyz_a[i_sol-1];
      status=fdevrul();
      if (status<0) 
         {
         sprintf(errbuf,"fdevrul (1)%%sur905*fanalys");
         return(varkon_erpush("SU2973",errbuf));
         }

      f_cur_p   = f_cur;
      dfdt_p    = dfdt;
/*!                                                                 */
/*    Recalculate function and derivative for - eps.                */
/*    Call of internal function fdevrul                             */
/*                                                                 !*/

      eps = - 0.00001;
      r_p = xyz_p.r;
      t_p = xyz_p.drdt;
      r_p.x_gm=xyz_p.r.x_gm+
           eps*xyz_p.drdt.x_gm+eps*eps/2*xyz_p.d2rdt2.x_gm;
      r_p.y_gm=xyz_p.r.y_gm+eps*xyz_p.drdt.y_gm+
           eps*eps/2*xyz_p.d2rdt2.y_gm;
      r_p.z_gm=xyz_p.r.z_gm+eps*xyz_p.drdt.z_gm+
           eps*eps/2*xyz_p.d2rdt2.z_gm;
      t_p.x_gm=xyz_p.drdt.x_gm+eps/2*xyz_p.d2rdt2.x_gm;
      t_p.y_gm=xyz_p.drdt.y_gm+eps/2*xyz_p.d2rdt2.y_gm;
      t_p.z_gm=xyz_p.drdt.z_gm+eps/2*xyz_p.d2rdt2.z_gm;
      xyz = xyz_a[i_sol-1];
      status=fdevrul();
      if (status<0) 
         {
         sprintf(errbuf,"fdevrul (1)%%sur905*fanalys");
         return(varkon_erpush("SU2973",errbuf));
         }


      f_cur_n   = f_cur;
      dfdt_n    = dfdt;


/*!                                                                 */
/*    Set the output code for the developable ruling:               */
/*    Eq. 1: Developable, normal surface for in both directions     */
/*    Eq. 2: Developable, normal surface in positive direction only */
/*    Eq. 3: Developable, normal surface in negative direction only */
/*                                                                 !*/

      d_code[i_sol-1] = I_UNDEF;
      if       (  f_cur_p >  0.0 && dfdt_p <  0.0 &&
                  f_cur_n >  0.0 && dfdt_n >  0.0    )
            d_code[i_sol-1] = 11;
      else if  (  f_cur_p <  0.0 && dfdt_p >  0.0 &&
                  f_cur_n <  0.0 && dfdt_n <  0.0    ) 
            d_code[i_sol-1] = 12;
      else if  (  f_cur_p >  0.0 && dfdt_p <  0.0 &&
                  f_cur_n <  0.0 && dfdt_n <  0.0    ) 
            d_code[i_sol-1] = 13;
      else if  (  f_cur_p <  0.0 && dfdt_p >  0.0 &&
                  f_cur_n >  0.0 && dfdt_n >  0.0    ) 
            d_code[i_sol-1] = 14;

      else if  (  f_cur_p <  0.0 && dfdt_p >  0.0 &&
                  f_cur_n >  0.0 && dfdt_n <  0.0    ) 
            d_code[i_sol-1] = 21;
      else if  (  f_cur_p <  0.0 && dfdt_p >  0.0 &&
                  f_cur_n <  0.0 && dfdt_n >  0.0    ) 
            d_code[i_sol-1] = 22;
      else if  (  f_cur_p >  0.0 && dfdt_p <  0.0 &&
                  f_cur_n <  0.0 && dfdt_n >  0.0    ) 
            d_code[i_sol-1] = 23;
      else if  (  f_cur_p >  0.0 && dfdt_p <  0.0 &&
                  f_cur_n >  0.0 && dfdt_n <  0.0    ) 
            d_code[i_sol-1] = 24;

      else if  (  f_cur_p <  0.0 && dfdt_p <  0.0 &&
                  f_cur_n >  0.0 && dfdt_n >  0.0    ) 
            d_code[i_sol-1] = 31;
      else if  (  f_cur_p >  0.0 && dfdt_p >  0.0 &&
                  f_cur_n <  0.0 && dfdt_n <  0.0    ) 
            d_code[i_sol-1] = 32;
      else if  (  f_cur_p <  0.0 && dfdt_p <  0.0 &&
                  f_cur_n <  0.0 && dfdt_n <  0.0    ) 
            d_code[i_sol-1] = 33;
      else if  (  f_cur_p >  0.0 && dfdt_p >  0.0 &&
                  f_cur_n >  0.0 && dfdt_n >  0.0    ) 
            d_code[i_sol-1] = 34;



      else if  (  f_cur_p <  0.0 && dfdt_p <  0.0 &&
                  f_cur_n >  0.0 && dfdt_n <  0.0    ) 
            d_code[i_sol-1] = 41;
      else if  (  f_cur_p >  0.0 && dfdt_p >  0.0 &&
                  f_cur_n <  0.0 && dfdt_n >  0.0    ) 
            d_code[i_sol-1] = 42;
      else if  (  f_cur_p <  0.0 && dfdt_p <  0.0 &&
                  f_cur_n <  0.0 && dfdt_n >  0.0    ) 
            d_code[i_sol-1] = 43;
      else if  (  f_cur_p >  0.0 && dfdt_p >  0.0 &&
                  f_cur_n >  0.0 && dfdt_n <  0.0    ) 
            d_code[i_sol-1] = 44;


      else
         {
            d_code[i_sol-1] = -1;
#ifdef  TILLSVIDARE
         sprintf(errbuf,"d_code%%sur905*fanalys");
         return(varkon_erpush("SU2993",errbuf));
#endif /*  TILLSVIDARE */
         }

#ifdef DEBUG
if ( dbglev(SURPAC) == 1  )
  {
  fprintf(dbgfil(SURPAC),
  "sur905*fanalys f_cur_p %12.2f dfdt_p %12.2f d_code %d\n",
   f_cur_p, dfdt_p,(int)d_code[i_sol-1]);
  }
if ( dbglev(SURPAC) == 1  )
  {
  fprintf(dbgfil(SURPAC),
  "sur905*fanalys f_cur_n %12.2f dfdt_n %12.2f \n",
   f_cur_n, dfdt_n);
  }
#endif

  if       ( d_code[i_sol-1] >= 11 && d_code[i_sol-1] <= 14 ) 
       d_code[i_sol-1] = 1;
  else if  ( d_code[i_sol-1] >= 21 && d_code[i_sol-1] <= 24 ) 
       d_code[i_sol-1] = 2;
  else if  ( d_code[i_sol-1] >= 31 && d_code[i_sol-1] <= 34 ) 
       d_code[i_sol-1] = 3;
  else if  ( d_code[i_sol-1] >= 41 && d_code[i_sol-1] <= 44 ) 
       d_code[i_sol-1] = 4;

#ifdef  KANSKE_FUNGERAR
/*    Derivative of f_cur for output point                          */
      dfdt_sol  = dfdt;
      f_cur_sol = f_cur;

/*!                                                                 */
/*    Let the solution point be r_p and t_p and xyz_p be the        */
/*    xyz derivatives for function fdevrul, which is called.        */
/*                                                                 !*/
      
      r_p = xyz_a[0].r;
      t_p = xyz_a[0].drdt;

      xyz = xyz_p;
      status=fdevrul();
      if (status<0) 
         {
         sprintf(errbuf,"fdevrul (2)%%sur905*fanalys");
         return(varkon_erpush("SU2973",errbuf));
         }

/*    Derivative when curves are switched                           */
      dfdt_switch = dfdt;

/*                                                                  */
/*    Set the output code for the developable ruling:               */
/*    Eq. 1: Derivatives both positive  <==> Normal developable     */
/*    Eq. 2: Derivatives both negative  <==> Normal developable     */
/*    Eq. 3: Derivative crv 1 positive  <==> Normal developable     */
/*    Eq. 4: Derivative crv 2 positive  <==> Normal developable     */
/*                                                                  */

      d_code[i_sol-1] = F_UNDEF;
      if       ( dfdt_sol >= 0.0 && dfdt_switch >= 0.0 ) 
            d_code[i_sol-1] = 1;
      else if  ( dfdt_sol <  0.0 && dfdt_switch <  0.0 ) 
            d_code[i_sol-1] = 2;
      else if  ( dfdt_sol >= 0.0 && dfdt_switch <  0.0 ) 
            d_code[i_sol-1] = 3;
      else if  ( dfdt_sol <  0.0 && dfdt_switch >= 0.0 ) 
            d_code[i_sol-1] = 4;
      else 
            d_code[i_sol-1] = -1;

#ifdef DEBUG
if ( dbglev(SURPAC) == 1  )
  {
  fprintf(dbgfil(SURPAC),
  "sur905*fanalys u %10.5f dfdt_sol %f dfdt_switch %f d_code %d\n",
   xyz_a[i_sol-1].t_global,dfdt_sol,dfdt_switch,(int)d_code[i_sol-1]);
  }
if ( dbglev(SURPAC) == 1  )
  {
  fprintf(dbgfil(SURPAC),
  "sur905*fanalys u %10.5f f_cur_sol %f f_cur_switch %f\n",
   xyz_a[i_sol-1].t_global,f_cur_sol,f_cur);
  }
#endif
#endif /*  KANSKE_FUNGERAR  */


     }  /* End loop i_sol= 1,2,..., numsol                          */



    return(SUCCED);


} /* End of function */

/*!****************************************************************!*/





Generated by  Doxygen 1.6.0   Back to index