/********************************************************************/ /* */ /* 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_sur_uvsegeval File: sur214.c */ /* ============================================================= */ /* */ /* Purpose */ /* ------- */ /* */ /* The function calculates the coordinates and derivatives */ /* for a point on a segment of a U,V (surface) curve */ /* */ /* Author: Gunnar Liden */ /* */ /* Revisions */ /* */ /* 1994-05-28 Originally written */ /* 1994-10-29 Debug */ /* 1994-11-12 Offset curve */ /* 1994-12-08 Codes for derivatives, debug */ /* 1995-02-25 u_t2_geod0, v_t2_geod0, Debug (DB pointer, ...) */ /* 1995-03-15 Documentation */ /* 1995-03-17 New DBread_one_patch(), J. Kjellander */ /* 1996-02-13 Numerical noise from evaluation of UV curve */ /* 1997-12-07 NURBS, temporary (first) implementation */ /* 1998-02-03 NURBS bug */ /* 1999-12-18 Free source code modifications */ /* */ /* */ /******************************************************************!*/ /* ------------- Short description of function -----------------*/ /* */ /*sdescr varkon_sur_uvsegeval Coord. and derivatives for UV pt */ /* */ /*------------------------------------------------------------- */ /*!---------------------- Theory -----------------------------------*/ /* Reference: Faux & Pratt p 110-111 and p 274 */ /* */ /*-----------------------------------------------------------------!*/ /* -------------- Function calls (internal) ------------------------*/ /* */ /* */ /*----------------------------------------------------------------- */ /* -- Static (common) variables for the functions in this file -----*/ /* */ /*----------------------------------------------------------------- */ /* -------------- Function calls (external) ------------------------*/ /* */ /* DBread_surface * Retrieve DB surface data */ /* varkon_sur_patadr * Patch adress for given U,V pt */ /* varkon_pat_eval * Evaluate patch */ /* varkon_seg_uveval * Evaluate UV segment */ /* varkon_seg_offset * Offset coord. & derivatives */ /* varkon_sur_normkappa * Normal curvature */ /* varkon_erpush * Error message to terminal */ /* */ /*----------------------------------------------------------------- */ /* ------------ Error messages and warnings ------------------------*/ /* */ /* SU2943 = Called function xxxx failed in varkon_sur_uvsegeval */ /* SU2993 = Severe program error in varkon_sur_uvsegeval (sur214). */ /* */ /*----------------------------------------------------------------- */ /*!****************** Function **************************************/ /* */ DBstatus varkon_sur_uvsegeval ( /*-------------- Argument declarations -----------------------------*/ /* */ /* In: */ DBCurve *p_cur, /* Curve (ptr) */ DBSeg *p_seg, /* Curve segment (ptr) */ EVALC *p_xyz_c ) /* Curve coordinates & derivatives (ptr) */ /* Out: */ /* Data to p_xyz_c. */ /*-----------------------------------------------------------------!*/ { /* Start of function */ /*!--------------- Internal variables ------------------------------*/ /* */ EVALS xyz_s; /* Surface coordinates & derivatives */ DBfloat dsdt; /* Length of tangent (dr/dt) */ DBfloat tan[3]; /* Normalized tangent= (dr/dt)/dsdt */ DBfloat nkappa; /* Normal curvature */ DBfloat geodesic; /* Geodesic curvature */ /* */ /*-----------------------------------------------------------------!*/ DBSurf sur; /* Surface for input U,V segment */ DBint icase; /* Evaluation case for varkon_pat_eval */ DBint nu,nv; /* Number of patches in surface */ DBint iu,iv; /* Patch number (address) */ DBfloat u_local; /* Local patch parameter U */ DBfloat v_local; /* Local patch parameter U */ DBVector drdt; /* First derivative with respect to t */ DBVector d2rdt2; /* Second derivative with respect to t */ DBfloat kappa; /* Curvature */ DBVector b_norm; /* Binormal */ DBVector p_norm; /* Principal normal */ bool offseg; /* Flag for offset segment */ DBfloat offset; /* Offset value */ DBVector n_p; /* Curve segment plane normal */ DBVector r; /* Point non-offset */ DBfloat u_t2_geod0; /* Second U derivative for geodesic= 0 */ DBfloat v_t2_geod0; /* Second V derivative for geodesic= 0 */ DBVector v1,v2,v3,v4; /* For geodesic= 0 (UV 2'nd derivatives) */ DBfloat dot_1,dot_2; /* Scalar products for geodesic = 0 */ DBfloat dot_3; /* Scalar product for geodesic = 0 */ DBfloat uv_t_len; /* Length u_t,v_t for geodesic = 0 */ DBint rcode; /* Derivative calculation code */ DBint status; /* Error code from called function */ char errbuf[80]; /* String for error message fctn erpush */ DBPatch *p_toppat; /* Pointer to topological patch */ DBPatch toppat; /* Topological patch for DBread_one_patch */ GMPATX geopat; /* Geometrical patch for DBread_one_patch */ #ifdef DEBUG DBfloat larserik; /* Geodesic curvature a la Lars-Erik */ DBfloat kvadrat; /* Geodesic curvature as .. */ #endif /*--------------end-of-declarations---------------------------------*/ /*!New-Page--------------------------------------------------------!*/ /*! !*/ /*! */ /* Algorithm */ /* ========= */ /* */ /* */ /* !*/ if ( p_xyz_c->evltyp & EVC_DR ) rcode = 1; #ifdef DEBUG if ( dbglev(SURPAC) == 1 && (p_xyz_c->evltyp & EVC_R) ) { fprintf(dbgfil(SURPAC), "sur214 Enter*varkon_sur_uvsegeval Case EVC_R Offset %f \n", p_seg->ofs); } else if ( dbglev(SURPAC) == 1 && (p_xyz_c->evltyp & EVC_DR) ) { fprintf(dbgfil(SURPAC), "sur214 Enter*varkon_sur_uvsegeval Case EVC_DR Offset %f \n", p_seg->ofs); } else { fprintf(dbgfil(SURPAC), "sur214 Enter*varkon_sur_uvsegeval Case= %d Offset %f \n", (int)p_xyz_c->evltyp , p_seg->ofs); fflush(dbgfil(SURPAC)); } #endif /*! */ /* 0. Check of input data and initializations */ /* ____________________________________________ */ /* */ /* Let offset flag be true and retrieve curve plane if the */ /* curve is in offset. */ /* !*/ /* Check that DB pointer to patch not is NULL (assume that DBSeg */ /* data is initialized with sur779) for Debug On */ #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 Patch pointer in DB p_seg->spek_gm %d\n", (int)p_seg->spek_gm ); fflush(dbgfil(SURPAC)); } if ( p_seg->spek_gm == DBNULL ) { sprintf(errbuf,"DB ptr is NULL%%sur214"); return(varkon_erpush("SU2993",errbuf)); } #endif /* Initialization of local variables */ iu = I_UNDEF; iv = I_UNDEF; nu = I_UNDEF; nv = I_UNDEF; u_local = F_UNDEF; v_local = F_UNDEF; if ( p_seg->ofs != 0.0 ) { offseg = TRUE; offset = p_seg->ofs; n_p.x_gm = p_cur->csy_cu.g31; n_p.y_gm = p_cur->csy_cu.g32; n_p.z_gm = p_cur->csy_cu.g33; } else { offseg = FALSE; offset = p_seg->ofs; n_p.x_gm = -0.123456789; n_p.y_gm = -0.123456789; n_p.z_gm = -0.123456789; } #ifdef DEBUG if ( dbglev(SURPAC) == 1 && offseg == TRUE) { fprintf(dbgfil(SURPAC), "sur214 Input segment is in offset= %f\n", offset); fprintf(dbgfil(SURPAC), "sur214 Curve plane normal n_p %f %f %f \n", n_p.x_gm, n_p.y_gm , n_p.z_gm); fflush(dbgfil(SURPAC)); } #endif /*! */ /*! */ /* 1. Surface U,V curve segment evaluation */ /* _______________________________________ */ /* */ /* Calculate coordinates u,v and derivatives du/dt,dvdt, ... */ /* for the given point on the input U,V segment. */ /* Call of varkon_seg_uveval (sur786). */ /* !*/ status=varkon_seg_uveval (p_cur,p_seg,p_xyz_c); if (status<0) { sprintf(errbuf,"sur768%%sur214"); return(varkon_erpush("SU2943",errbuf)); } /*! */ /* 2. Current patch address and local parameter point. */ /* ___________________________________________________ */ /* */ /* Retrieve surface header data (no of patches) from DB. */ /* Call of DBread_surface. */ /* !*/ status= DBread_surface (&sur,p_seg->spek_gm); if (status<0) { sprintf(errbuf,"DBread_surface%%varkon_sur_uvsegeval"); return(varkon_erpush("SU2943",errbuf)); } nu = sur.nu_su; nv = sur.nv_su; #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 Surface data retrieved from DB with DBread_surface nu %d nv %d\n", (int)nu , (int)nv ); fflush(dbgfil(SURPAC)); } #endif /*! */ /* Patch address iu,iv and local patch U,V point. */ /* Eliminate numerical noise from evaluation of UV curve */ /* Call of varkon_sur_patadr (sur211). */ /* !*/ #ifdef DEBUG if ( sur.typ_su != NURB_SUR ) { if ( dbglev(SURPAC) == 1 && p_xyz_c->u < 1.0 ) { fprintf(dbgfil(SURPAC), "sur214 Evaluated u= %25.15f Modified to 1.0 (Diff. %25.15f)\n", p_xyz_c->u, p_xyz_c->u - 1.0 ); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 1 && p_xyz_c->v < 1.0 ) { fprintf(dbgfil(SURPAC), "sur214 Evaluated v= %25.15f Modified to 1.0 (Diff. %25.15f)\n", p_xyz_c->v, p_xyz_c->v - 1.0 ); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 1 && p_xyz_c->u > (DBfloat)nu + 1.0 ) { fprintf(dbgfil(SURPAC), "sur214 Evaluated u= %25.15f Modified to 1.0 (Diff. %25.15f)\n", p_xyz_c->u, p_xyz_c->u - (DBfloat)nu - 1.0 ); fflush(dbgfil(SURPAC)); } if ( dbglev(SURPAC) == 1 && p_xyz_c->v > (DBfloat)nv + 1.0 ) { fprintf(dbgfil(SURPAC), "sur214 Evaluated u= %25.15f Modified to 1.0 (Diff. %25.15f)\n", p_xyz_c->v, p_xyz_c->v - (DBfloat)nv - 1.0 ); fflush(dbgfil(SURPAC)); } } /* End Not NURBS */ #endif if ( sur.typ_su != NURB_SUR ) { if ( p_xyz_c->u < 1.0 ) p_xyz_c->u = 1.0; if ( p_xyz_c->v < 1.0 ) p_xyz_c->v = 1.0; if ( p_xyz_c->u > (DBfloat)nu + 1.0 ) p_xyz_c->u = (DBfloat)nu+1.0; if ( p_xyz_c->v > (DBfloat)nv + 1.0 ) p_xyz_c->v = (DBfloat)nv+1.0; status=varkon_sur_patadr (p_xyz_c->u,p_xyz_c->v,nu,nv,&iu,&iv,&u_local,&v_local); if (status<0) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 sur211 failed for u %25.15f v %25.15f \n", p_xyz_c->u, p_xyz_c->v ); fflush(dbgfil(SURPAC)); } #endif sprintf(errbuf,"varkon_sur_patadr%%varkon_sur_uvsegeval"); return(varkon_erpush("SU2943",errbuf)); } } /* End Not NURBS */ /*! */ /* 3. Surface coordinates and derivatives */ /* ______________________________________ */ /* */ /* Get C pointer to patch (retrieve patch data). */ /* Call of DBread_one_patch. */ /* !*/ /* A NURBS surface consists of only one patch and cannot be */ /* partially be retrieved. It is for all other surface types */ /* only necessary to retrieve one patch */ switch ( sur.typ_su ) { case NURB_SUR: DBread_patches( &sur,&p_toppat); u_local = p_xyz_c->u; v_local = p_xyz_c->v; #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 NURBS First topology and the geometry patch retrieved\n"); fprintf(dbgfil(SURPAC), "sur214 Parameter point from sur786 EVALC calculation %f %f \n", u_local, v_local); fflush(dbgfil(SURPAC)); } #endif break; default: toppat.spek_c = (char *)&geopat; if ( (status=DBread_one_patch ( &sur,&toppat,(short)iu,(short)iv)) < 0 ) return(status); p_toppat = &toppat; break; } /*! */ /* Calculate coordinates and derivatives in the patch. */ /* Call of varkon_pat_eval (sur220). */ /* !*/ /* TODO: */ /* If the curve segment is an (iso) V- or U-segment can the */ /* the calculation time be reduced if icase is set to 1 or 2 */ #ifdef TODO_INVESTIGATE_FIRST if ( p_seg->c1y==0.0 && p_seg->c2y==0.0 && p_seg->c3y==0.0 && rcode==1 ) icase = 1; else if ( p_seg->c1x==0.0 && p_seg->c2x==0.0 && p_seg->c3y==0.0 && rcode==1 ) icase = 2; else icase = 3; #endif icase = 3; if (offseg == TRUE ) icase = 3; /* All derivatives for offset */ #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 Parameter for varkon_pat_eval u_local %f v_local %f \n", u_local, v_local); fflush(dbgfil(SURPAC)); } #endif status= varkon_pat_eval (&sur,p_toppat ,icase,u_local,v_local, &xyz_s); if (status<0) { /* Deallocate NURBS */ if ( sur.typ_su == NURB_SUR ) DBfree_patches( &sur,p_toppat); sprintf(errbuf,"varkon_pat_eval%%varkon_sur_uvsegeval"); return(varkon_erpush("SU2943",errbuf)); } /* Memory must be freed for a NURBS surface */ if ( sur.typ_su == NURB_SUR ) DBfree_patches( &sur,p_toppat); #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 p_xyz_c->u_t %f v_t %f u_t2 %f v_t2 %f\n", p_xyz_c->u_t, p_xyz_c->v_t ,p_xyz_c->u_t2,p_xyz_c->v_t2); } if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 xyz_s= %f %f %f \n", xyz_s.r_x , xyz_s.r_y ,xyz_s.r_z ); fflush(dbgfil(SURPAC)); } #endif /*! */ /* Coordinates from xyz_s. Goto nomore if evltyp= EVC_R. */ /* Goto nomore if evltyp= EVC_R and not an offset curve */ /* !*/ rcode = -1; if ( p_xyz_c->evltyp & EVC_R ) rcode = 0; if ( p_xyz_c->evltyp & EVC_DR ) rcode = 1; if ( p_xyz_c->evltyp & EVC_D2R ) rcode = 2; if ( (p_xyz_c->evltyp & EVC_PN ) || (p_xyz_c->evltyp & EVC_BN ) || (p_xyz_c->evltyp & EVC_KAP) ) rcode = 3; if ( rcode < 0 ) { sprintf(errbuf,"rcode%%varkon_sur_uvsegeval"); return(varkon_erpush("SU2993",errbuf)); } #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 rcode %d offseg %d (1)\n", (int)rcode, (int)offseg ); fflush(dbgfil(SURPAC)); } #endif p_xyz_c->r.x_gm= xyz_s.r_x; p_xyz_c->r.y_gm= xyz_s.r_y; p_xyz_c->r.z_gm= xyz_s.r_z; { if ( rcode == 0 && offseg == FALSE ) goto nomore; } r.x_gm= xyz_s.r_x; /* For offset calculation sur217 */ r.y_gm= xyz_s.r_y; r.z_gm= xyz_s.r_z; /*! */ /* First derivatives to output variable drdt. */ /* Goto nomore if evltyp= EVC_DR and not an offset curve */ /* !*/ drdt.x_gm= xyz_s.u_x*p_xyz_c->u_t + xyz_s.v_x*p_xyz_c->v_t; drdt.y_gm= xyz_s.u_y*p_xyz_c->u_t + xyz_s.v_y*p_xyz_c->v_t; drdt.z_gm= xyz_s.u_z*p_xyz_c->u_t + xyz_s.v_z*p_xyz_c->v_t; #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur214 rcode %d offseg %d (2)\n", (int)rcode, (int)offseg ); } #endif p_xyz_c->drdt.x_gm= drdt.x_gm; p_xyz_c->drdt.y_gm= drdt.y_gm; p_xyz_c->drdt.z_gm= drdt.z_gm; if ( rcode == 1 && offseg == FALSE ) { #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur214*Exit*varkon_sur_uvsegeval Output tangent %f %f %f\n", p_xyz_c->drdt.x_gm,p_xyz_c->drdt.y_gm,p_xyz_c->drdt.z_gm); } #endif goto nomore; } /*! */ /* Second derivatives to output variable d2rdt2 */ /* Goto nomore if evltyp= EVC_D2R and not an offset curve */ /* !*/ d2rdt2.x_gm= xyz_s.u2_x*p_xyz_c->u_t*p_xyz_c->u_t + 2.0* xyz_s.uv_x*p_xyz_c->u_t*p_xyz_c->v_t + xyz_s.v2_x*p_xyz_c->v_t*p_xyz_c->v_t + xyz_s.u_x *p_xyz_c->u_t2 + xyz_s.v_x *p_xyz_c->v_t2; d2rdt2.y_gm= xyz_s.u2_y*p_xyz_c->u_t*p_xyz_c->u_t + 2.0* xyz_s.uv_y*p_xyz_c->u_t*p_xyz_c->v_t + xyz_s.v2_y*p_xyz_c->v_t*p_xyz_c->v_t + xyz_s.u_y *p_xyz_c->u_t2 + xyz_s.v_y *p_xyz_c->v_t2; d2rdt2.z_gm= xyz_s.u2_z*p_xyz_c->u_t*p_xyz_c->u_t + 2.0* xyz_s.uv_z*p_xyz_c->u_t*p_xyz_c->v_t + xyz_s.v2_z*p_xyz_c->v_t*p_xyz_c->v_t + xyz_s.u_z *p_xyz_c->u_t2 + xyz_s.v_z *p_xyz_c->v_t2; p_xyz_c->d2rdt2.x_gm= d2rdt2.x_gm; p_xyz_c->d2rdt2.y_gm= d2rdt2.y_gm; p_xyz_c->d2rdt2.z_gm= d2rdt2.z_gm; if ( rcode == 2 && offseg == FALSE ) { #ifdef DEBUG if ( dbglev(SURPAC) == 2 ) { fprintf(dbgfil(SURPAC), "sur214*Exit*varkon_sur_uvsegeval Output d2rdt2 %f %f %f\n", p_xyz_c->d2rdt2.x_gm,p_xyz_c->d2rdt2.y_gm,p_xyz_c->d2rdt2.z_gm); } #endif goto nomore; } /*! */ /* Calculate offset coordinates and derivatives if the curve */ /* segment is in offset. */ /* Call of varkon_seg_offset (sur217). */ /* !*/ if (offseg == TRUE ) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 Before sur217 r %f %f %f\n", r.x_gm,r.y_gm,r.z_gm); fprintf(dbgfil(SURPAC), "sur214 Before sur217 drdt %f %f %f\n", drdt.x_gm,drdt.y_gm,drdt.z_gm); fprintf(dbgfil(SURPAC), "sur214 Before sur217 d2rdt2 %f %f %f\n", d2rdt2.x_gm,d2rdt2.y_gm,d2rdt2.z_gm); fflush(dbgfil(SURPAC)); } #endif status= varkon_seg_offset (r,drdt,d2rdt2,offset,p_xyz_c->evltyp,n_p, &r,&drdt,&d2rdt2); #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 After sur217 r %f %f %f\n", r.x_gm,r.y_gm,r.z_gm); fprintf(dbgfil(SURPAC), "sur214 After sur217 drdt %f %f %f\n", drdt.x_gm,drdt.y_gm,drdt.z_gm); fprintf(dbgfil(SURPAC), "sur214 After sur217 d2rdt2 %f %f %f\n", d2rdt2.x_gm,d2rdt2.y_gm,d2rdt2.z_gm); fprintf(dbgfil(SURPAC), "sur214 After sur217 status %d\n", (int)status ); fflush(dbgfil(SURPAC)); } #endif if (status<0) { sprintf(errbuf,"varkon_seg_offset%%varkon_sur_uvsegeval"); return(varkon_erpush("SU2943",errbuf)); } /* This must be wrong ????? TODO { */ p_xyz_c->r.x_gm= r.x_gm; p_xyz_c->r.y_gm= r.y_gm; p_xyz_c->r.z_gm= r.z_gm; if ( rcode == 0 ) { goto nomore; } p_xyz_c->drdt.x_gm= drdt.x_gm; p_xyz_c->drdt.y_gm= drdt.y_gm; p_xyz_c->drdt.z_gm= drdt.z_gm; if ( rcode == 1 ) { #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214*Exit*varkon_sur_uvsegeval Output tangent %f %f %f\n", p_xyz_c->drdt.x_gm,p_xyz_c->drdt.y_gm,p_xyz_c->drdt.z_gm); fflush(dbgfil(SURPAC)); } #endif goto nomore; } { p_xyz_c->d2rdt2.x_gm= d2rdt2.x_gm; p_xyz_c->d2rdt2.y_gm= d2rdt2.y_gm; p_xyz_c->d2rdt2.z_gm= d2rdt2.z_gm; if ( rcode == 2 ) #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214*Exit*varkon_sur_uvsegeval Output tangent %f %f %f\n", p_xyz_c->d2rdt2.x_gm,p_xyz_c->d2rdt2.y_gm,p_xyz_c->d2rdt2.z_gm); fflush(dbgfil(SURPAC)); } #endif goto nomore; } } /* End offseg == TRUE */ /*!New-Page--------------------------------------------------------!*/ /*! !*/ /*! !*/ /*! !*/ /*! */ /* 4. Calculate kappa and binormal */ /* _______________________________ */ /* !*/ /*! */ /* The length dsdt of the tangent (dr/dt) */ /* !*/ dsdt= SQRT(drdt.x_gm*drdt.x_gm + drdt.y_gm*drdt.y_gm + drdt.z_gm*drdt.z_gm); if ( dsdt < 1e-10 ) dsdt= 1e-10; #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 drdt %10.7f %10.7f %10.7f dsdt %f\n", drdt.x_gm,drdt.y_gm,drdt.z_gm ,dsdt); fflush(dbgfil(SURPAC)); } #endif /*! */ /* First kappa*binormal. Then kappa which is the length of */ /* vector kappa*binormal and lastly the binormal (divide */ /* vector kappa*binormal with kappa). */ /* !*/ b_norm.x_gm= (drdt.y_gm*d2rdt2.z_gm-drdt.z_gm*d2rdt2.y_gm)/dsdt/dsdt/dsdt; b_norm.y_gm= (drdt.z_gm*d2rdt2.x_gm-drdt.x_gm*d2rdt2.z_gm)/dsdt/dsdt/dsdt; b_norm.z_gm= (drdt.x_gm*d2rdt2.y_gm-drdt.y_gm*d2rdt2.x_gm)/dsdt/dsdt/dsdt; kappa = SQRT(b_norm.x_gm*b_norm.x_gm+ b_norm.y_gm*b_norm.y_gm+ b_norm.z_gm*b_norm.z_gm); p_xyz_c->kappa = kappa; if ( kappa > 1e-10 ) { b_norm.x_gm= b_norm.x_gm/kappa; b_norm.y_gm= b_norm.y_gm/kappa; b_norm.z_gm= b_norm.z_gm/kappa; } else { b_norm.x_gm= 0.0; b_norm.y_gm= 0.0; b_norm.z_gm= 0.0; } p_xyz_c->b_norm= b_norm; /*! */ /* 5. Calculate the principal normal */ /* _________________________________ */ /* !*/ /*! */ /* The normalized tangent */ /* !*/ tan[0] = drdt.x_gm/dsdt; tan[1] = drdt.y_gm/dsdt; tan[2] = drdt.z_gm/dsdt; /*! */ /* The principal normal= binormal X tangent */ /* !*/ p_norm.x_gm= b_norm.y_gm*tan[2]-b_norm.z_gm*tan[1]; p_norm.y_gm= b_norm.z_gm*tan[0]-b_norm.x_gm*tan[2]; p_norm.z_gm= b_norm.x_gm*tan[1]-b_norm.y_gm*tan[0]; p_xyz_c->p_norm= p_norm; /*! */ /* 6. Calculate the normal curvature */ /* _________________________________ */ /* */ /* Calculate the normal curvature */ /* Call of varkon_sur_normkappa (sur952). */ /* !*/ varkon_sur_normkappa (p_xyz_c->u_t,p_xyz_c->v_t,&xyz_s,&nkappa); /* No errors from this function */ p_xyz_c->nkappa = nkappa; /*! */ /* 7. Geodesic curvature */ /* _____________________ */ /* */ /* Calculate the geodesic (tangential) curvature */ /* !*/ geodesic=((drdt.y_gm*d2rdt2.z_gm-drdt.z_gm*d2rdt2.y_gm)*xyz_s.n_x+ (drdt.z_gm*d2rdt2.x_gm-drdt.x_gm*d2rdt2.z_gm)*xyz_s.n_y+ (drdt.x_gm*d2rdt2.y_gm-drdt.y_gm*d2rdt2.x_gm)*xyz_s.n_z) /dsdt/dsdt/dsdt; p_xyz_c->geodesic = geodesic; /*! */ /* Calculate the U,V second derivatives for geodesic curvature = 0 */ /* Can be used for calculation of geodesic curves or planar */ /* intersect curves (tangential curvature is zero for such curve) */ /* !*/ v1.x_gm= xyz_s.u2_x*p_xyz_c->u_t*p_xyz_c->u_t + 2.0* xyz_s.uv_x*p_xyz_c->u_t*p_xyz_c->v_t + xyz_s.v2_x*p_xyz_c->v_t*p_xyz_c->v_t; v1.y_gm= xyz_s.u2_y*p_xyz_c->u_t*p_xyz_c->u_t + 2.0* xyz_s.uv_y*p_xyz_c->u_t*p_xyz_c->v_t + xyz_s.v2_y*p_xyz_c->v_t*p_xyz_c->v_t; v1.z_gm= xyz_s.u2_z*p_xyz_c->u_t*p_xyz_c->u_t + 2.0* xyz_s.uv_z*p_xyz_c->u_t*p_xyz_c->v_t + xyz_s.v2_z*p_xyz_c->v_t*p_xyz_c->v_t; v4.x_gm = tan[1]* xyz_s.n_z-tan[2]* xyz_s.n_y; v4.y_gm = tan[2]* xyz_s.n_x-tan[0]* xyz_s.n_z; v4.z_gm = tan[0]* xyz_s.n_y-tan[1]* xyz_s.n_x; dot_1 = v1.x_gm*v4.x_gm + v1.y_gm*v4.y_gm + v1.z_gm*v4.z_gm; v2.x_gm = xyz_s.u_x; v2.y_gm = xyz_s.u_y; v2.z_gm = xyz_s.u_z; v3.x_gm = xyz_s.v_x; v3.y_gm = xyz_s.v_y; v3.z_gm = xyz_s.v_z; dot_2 = v2.x_gm*v4.x_gm + v2.y_gm*v4.y_gm + v2.z_gm*v4.z_gm; dot_3 = v3.x_gm*v4.x_gm + v3.y_gm*v4.y_gm + v3.z_gm*v4.z_gm; uv_t_len = SQRT(p_xyz_c->u_t*p_xyz_c->u_t + p_xyz_c->v_t*p_xyz_c->v_t); if ( fabs(dot_2) > fabs(dot_3) ) { v_t2_geod0 = uv_t_len*0.5; if ( fabs(dot_2) > 0.00000000001 ) { u_t2_geod0 = (- dot_1 - v_t2_geod0*dot_3)/dot_2; } else { sprintf(errbuf,"dot_2= 0 %%varkon_sur_uvsegeval"); return(varkon_erpush("SU2993",errbuf)); } } else { u_t2_geod0 = uv_t_len*0.5; if ( fabs(dot_3) > 0.00000000001 ) { v_t2_geod0 = (- dot_1 - u_t2_geod0*dot_2)/dot_3; } else { sprintf(errbuf,"dot_3= 0 %%varkon_sur_uvsegeval"); return(varkon_erpush("SU2993",errbuf)); } } p_xyz_c->u_t2_geod0 = u_t2_geod0; p_xyz_c->v_t2_geod0 = v_t2_geod0; #ifdef DEBUG v1.x_gm= xyz_s.u2_x*p_xyz_c->u_t*p_xyz_c->u_t + 2.0* xyz_s.uv_x*p_xyz_c->u_t*p_xyz_c->v_t + xyz_s.v2_x*p_xyz_c->v_t*p_xyz_c->v_t; v1.y_gm= xyz_s.u2_y*p_xyz_c->u_t*p_xyz_c->u_t + 2.0* xyz_s.uv_y*p_xyz_c->u_t*p_xyz_c->v_t + xyz_s.v2_y*p_xyz_c->v_t*p_xyz_c->v_t; v1.z_gm= xyz_s.u2_z*p_xyz_c->u_t*p_xyz_c->u_t + 2.0* xyz_s.uv_z*p_xyz_c->u_t*p_xyz_c->v_t + xyz_s.v2_z*p_xyz_c->v_t*p_xyz_c->v_t; v2.x_gm= xyz_s.u_x*p_xyz_c->u_t2; v2.y_gm= xyz_s.u_y*p_xyz_c->u_t2; v2.z_gm= xyz_s.u_z*p_xyz_c->u_t2; v3.x_gm= xyz_s.v_x*p_xyz_c->v_t2; v3.y_gm= xyz_s.v_y*p_xyz_c->v_t2; v3.z_gm= xyz_s.v_z*p_xyz_c->v_t2; v2.x_gm= xyz_s.u_x*p_xyz_c->u_t2_geod0; /* Test geodesic = 0 */ v2.y_gm= xyz_s.u_y*p_xyz_c->u_t2_geod0; v2.z_gm= xyz_s.u_z*p_xyz_c->u_t2_geod0; v3.x_gm= xyz_s.v_x*p_xyz_c->v_t2_geod0; v3.y_gm= xyz_s.v_y*p_xyz_c->v_t2_geod0; v3.z_gm= xyz_s.v_z*p_xyz_c->v_t2_geod0; v4.x_gm = tan[1]* xyz_s.n_z-tan[2]* xyz_s.n_y; v4.y_gm = tan[2]* xyz_s.n_x-tan[0]* xyz_s.n_z; v4.z_gm = tan[0]* xyz_s.n_y-tan[1]* xyz_s.n_x; larserik= v1.x_gm*v4.x_gm + v1.y_gm*v4.y_gm + v1.z_gm*v4.z_gm + v2.x_gm*v4.x_gm + v2.y_gm*v4.y_gm + v2.z_gm*v4.z_gm + v3.x_gm*v4.x_gm + v3.y_gm*v4.y_gm + v3.z_gm*v4.z_gm; larserik= larserik/dsdt/dsdt; kvadrat = kappa*kappa - nkappa*nkappa; if ( kvadrat > 0.0 ) kvadrat = SQRT(kvadrat); else kvadrat = -1.23456789; if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214 geodesic curvature %20.16f larserik %20.16f\n", geodesic, larserik ); fflush(dbgfil(SURPAC)); fprintf(dbgfil(SURPAC), "sur214 kappa %20.16f nkappa %20.16f\n", kappa , nkappa ); fflush(dbgfil(SURPAC)); fprintf(dbgfil(SURPAC), "sur214 geodesic kvadrat %20.16f \n", kvadrat ); fflush(dbgfil(SURPAC)); } #endif nomore: /*! Label nomore: End of calculation for case evltyp= .... !*/ #ifdef DEBUG if ( dbglev(SURPAC) == 1 ) { fprintf(dbgfil(SURPAC), "sur214*Exit*varkon_sur_uvsegeval Output pt %f %f %f\n", p_xyz_c->r.x_gm,p_xyz_c->r.y_gm,p_xyz_c->r.z_gm); fflush(dbgfil(SURPAC)); } #endif return(SUCCED); } /* End of function */ /********************************************************************/