ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/UseTheForce/Shapes_FF.cpp
(Generate patch)

Comparing trunk/OOPSE-2.0/src/UseTheForce/Shapes_FF.cpp (file contents):
Revision 1646 by chrisfen, Tue Oct 26 18:02:46 2004 UTC vs.
Revision 1688 by chrisfen, Fri Oct 29 22:28:12 2004 UTC

# Line 11 | Line 11
11   #include <stdio.h>
12   #include <string.h>
13   #include <map>
14 + #include <cmath>
15   #include <iostream>
16  
16 – using namespace std;
17 – using namespace oopse;
18 –
17   #ifdef IS_MPI
18   #include <mpi.h>
19   #endif //is_mpi
# Line 23 | Line 21 | using namespace oopse;
21   #include "UseTheForce/ForceFields.hpp"
22   #include "primitives/SRI.hpp"
23   #include "utils/simError.h"
24 + #include "utils/StringUtils.hpp"
25   #include "io/basic_ifstrstream.hpp"
26   #include "math/RealSphericalHarmonic.hpp"
27   #include "math/SquareMatrix3.hpp"
28 <
28 > #include "types/ShapeAtomType.hpp"
29   #include "UseTheForce/DarkSide/atype_interface.h"
30   #include "UseTheForce/DarkSide/shapes_interface.h"
31  
# Line 34 | Line 33 | Shapes_FF::Shapes_FF() {
33   #include "UseTheForce/mpiForceField.h"
34   #endif // is_mpi
35  
36 < Shapes_FF::Shapes_FF() {
37 <  Shapes_FF("");
39 < }
36 > using namespace std;
37 > using namespace oopse;
38  
41 – Shapes_FF::Shapes_FF(char* the_variant){
42 –  ffPath_env = "FORCE_PARAM_PATH";
43 –
44 – }
45 –
39   Shapes_FF::~Shapes_FF(){
40  
41   #ifdef IS_MPI
42    if( worldRank == 0 ){
43   #endif // is_mpi
44      
45 <    fclose( frcFile );
45 >    forceFile.close();
46      
47   #ifdef IS_MPI
48    }
# Line 60 | Line 53 | void Shapes_FF::calcRcut( void ){
53   void Shapes_FF::calcRcut( void ){
54    
55   #ifdef IS_MPI
56 <  double tempShapesRcut = shapesRcut;
56 >  double tempShapesRcut = bigContact;
57    MPI_Allreduce( &tempShapesRcut, &shapesRcut, 1, MPI_DOUBLE, MPI_MAX,
58                   MPI_COMM_WORLD);
59 + #else
60 +  shapesRcut = bigContact;
61   #endif  //is_mpi
62    entry_plug->setDefaultRcut(shapesRcut);
63   }
# Line 73 | Line 68 | void Shapes_FF::readForceFile( void ){
68   }
69  
70  
71 < void Shapes_FF::readForceFile( void ){
71 > void Shapes_FF::readParams( void ){
72  
73    char readLine[1024];
74 <  char fileName[200];
75 <  char temp[200];
76 <  char* ffPath;
77 <  char *shapeFileName;
74 >
75 >  string fileName;
76 >  string shapeFileName;
77 >  string tempString;
78 >
79    char *nameToken;
80    char *delim = " ,;\t\n";
81    int nTokens, i;
# Line 87 | Line 83 | void Shapes_FF::readForceFile( void ){
83    int nRange = 0;
84    int nStrength = 0;
85    int myATID;
86 +  int isError;
87    string nameString;
88 <  ShapeType* st;
89 <  map<string, AtomType*> atomTypeMap;
88 >  AtomType* at;
89 >  DirectionalAtomType* dat;
90 >  ShapeAtomType* st;
91 >
92    map<string, AtomType*>::iterator iter;
93  
94    // vectors for shape transfer to fortran
95 <  vector<RealSphericalHarmonic> tempSHVector;
95 >  vector<RealSphericalHarmonic*> tempSHVector;
96    vector<int> contactL;
97    vector<int> contactM;
98    vector<int> contactFunc;
# Line 106 | Line 105 | void Shapes_FF::readForceFile( void ){
105    vector<int> strengthM;
106    vector<int> strengthFunc;
107    vector<double> strengthCoeff;
109 –
110 –  // build a basic file reader
111 –  ifstrsteam frcReader;
108    
109    // generate the force file name  
110 <  strcpy( fileName, "Shapes.frc" );
110 >  fileName = "Shapes.frc";
111    
112    // attempt to open the file in the current directory first.
113 <  frcReader.open( fileName );
113 >  forceFile.open( fileName.c_str() );
114    
115 <  if( frcReader == NULL ){
115 >  if( forceFile == NULL ){
116      
117 <    // next see if the force path enviorment variable is set
117 >    tempString = ffPath;
118 >    tempString += "/";
119 >    tempString += fileName;
120 >    fileName = tempString;
121      
122 <    ffPath = getenv( ffPath_env );
124 <    if( ffPath == NULL ) {
125 <      STR_DEFINE(ffPath, FRC_PATH );
126 <    }
127 <        
128 <    strcpy( temp, ffPath );
129 <    strcat( temp, "/" );
130 <    strcat( temp, fileName );
131 <    strcpy( fileName, temp );
122 >    forceFile.open( fileName.c_str() );
123      
124 <    frcReader.open( fileName );
134 <    
135 <    if( frcFile == NULL ){
124 >    if( forceFile == NULL ){
125        
126        sprintf( painCave.errMsg,
127                 "Error opening the force field parameter file:\n"
128                 "\t%s\n"
129                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
130                 "variable?\n",
131 <               fileName );
131 >               fileName.c_str() );
132        painCave.severity = OOPSE_ERROR;
133        painCave.isFatal = 1;
134        simError();
# Line 148 | Line 137 | void Shapes_FF::readForceFile( void ){
137    
138    // read in the shape types.
139    
140 <  findBegin( "ShapeTypes" );
140 >  findBegin( forceFile, "ShapeTypes" );
141    
142 <  while( !frcReader.eof() ){
143 <    frcReader.getline( readLine, sizeof(readLine) );
142 >  while( !forceFile.eof() ){
143 >    forceFile.getline( readLine, sizeof(readLine) );
144  
145      // toss comment lines
146      if( readLine[0] != '!' && readLine[0] != '#' ){
147        
148        if (isEndLine(readLine)) break;
149        
150 <      nTokens = count_tokens(readLine, " ,;\t");
150 >      nTokens = countTokens(readLine, " ,;\t");
151        if (nTokens != 0) {
152          if (nTokens < 2) {
153            sprintf( painCave.errMsg,
154 <                   "Not enough data on a ShapeTypes line in file: %s\n"
155 <                   fileName );
154 >                   "Not enough data on a ShapeTypes line in file: %s\n",
155 >                   fileName.c_str() );
156            painCave.severity = OOPSE_ERROR;
157            painCave.isFatal = 1;
158            simError();
# Line 180 | Line 169 | void Shapes_FF::readForceFile( void ){
169          if (iter == atomTypeMap.end()) {
170            // no, it doesn't, so we may proceed:
171            
172 <          st = new ShapeType();
172 >          st = new ShapeAtomType();
173  
174            st->setName(nameString);
175 <          myATID = atomTypeMap.size();
175 >          myATID = atomTypeMap.size() + 1;
176            st->setIdent(myATID);
177            parseShapeFile(shapeFileName, st);
178            st->complete();
# Line 194 | Line 183 | void Shapes_FF::readForceFile( void ){
183            // declared in a previous block, and we just need to add
184            // the shape-specific information for this AtomType:
185  
186 <          st = (ShapeType*)(iter->second);
186 >          st = (ShapeAtomType*)(iter->second);
187            parseShapeFile(shapeFileName, st);
188          }
189        }
# Line 213 | Line 202 | void Shapes_FF::readForceFile( void ){
202    // pack up and send the necessary info to fortran
203    for (iter = atomTypeMap.begin(); iter != atomTypeMap.end(); ++iter){
204  
205 <    at = (AtomType*)(iter.second);
205 >    at = (AtomType*)(iter->second);
206  
207      if (at->isDirectional()) {
208  
# Line 232 | Line 221 | void Shapes_FF::readForceFile( void ){
221          
222          nContact = tempSHVector.size();
223          for (i=0; i<nContact; i++){
224 <          contactL.push_back(tempSHVector[i].getL());
225 <          contactM.push_back(tempSHVector[i].getM());
226 <          contactFunc.push_back(tempSHVector[i].getFunctionType());
227 <          contactCoeff.push_back(tempSHVector[i].getCoefficient());
224 >          contactL.push_back(tempSHVector[i]->getL());
225 >          contactM.push_back(tempSHVector[i]->getM());
226 >          contactFunc.push_back(tempSHVector[i]->getFunctionType());
227 >          contactCoeff.push_back(tempSHVector[i]->getCoefficient());
228          }
229          
230          rangeL.clear();
# Line 247 | Line 236 | void Shapes_FF::readForceFile( void ){
236          
237          nRange = tempSHVector.size();
238          for (i=0; i<nRange; i++){
239 <          rangeL.push_back(tempSHVector[i].getL());
240 <          rangeM.push_back(tempSHVector[i].getM());
241 <          rangeFunc.push_back(tempSHVector[i].getFunctionType());
242 <          rangeCoeff.push_back(tempSHVector[i].getCoefficient());
239 >          rangeL.push_back(tempSHVector[i]->getL());
240 >          rangeM.push_back(tempSHVector[i]->getM());
241 >          rangeFunc.push_back(tempSHVector[i]->getFunctionType());
242 >          rangeCoeff.push_back(tempSHVector[i]->getCoefficient());
243          }
244          
245          strengthL.clear();
# Line 262 | Line 251 | void Shapes_FF::readForceFile( void ){
251          
252          nStrength = tempSHVector.size();
253          for (i=0; i<nStrength; i++){
254 <          strengthL.push_back(tempSHVector[i].getL());
255 <          strengthM.push_back(tempSHVector[i].getM());
256 <          strengthFunc.push_back(tempSHVector[i].getFunctionType());
257 <          strengthCoeff.push_back(tempSHVector[i].getCoefficient());
254 >          strengthL.push_back(tempSHVector[i]->getL());
255 >          strengthM.push_back(tempSHVector[i]->getM());
256 >          strengthFunc.push_back(tempSHVector[i]->getFunctionType());
257 >          strengthCoeff.push_back(tempSHVector[i]->getCoefficient());
258          }
259          
260          isError = 0;
261          myATID = at->getIdent();
262 <        makeShape( &nContact, &contactL, &contactM, &contactFunc,
263 <                   &contactCoeff,
264 <                   &nRange, &rangeL, &rangeM, &rangeFunc, &rangeCoeff,
265 <                   &nStrength, &strengthL, &strengthM,
266 <                   &strengthFunc, &strengthCoeff,
262 >
263 >        makeShape( &nContact, &contactL[0], &contactM[0], &contactFunc[0],
264 >                   &contactCoeff[0],
265 >                   &nRange, &rangeL[0], &rangeM[0], &rangeFunc[0],
266 >                   &rangeCoeff[0],
267 >                   &nStrength, &strengthL[0], &strengthM[0],
268 >                   &strengthFunc[0], &strengthCoeff[0],
269                     &myATID,
270                     &isError);
271 +
272          if( isError ){
273            sprintf( painCave.errMsg,
274                     "Error initializing the \"%s\" shape in fortran\n",
275 <                   marker->first );
275 >                   (iter->first).c_str() );
276            painCave.isFatal = 1;
277            simError();
278          }
# Line 288 | Line 280 | void Shapes_FF::readForceFile( void ){
280      }
281    }
282    
283 +  isError = 0;
284 +  completeShapeFF(&isError);
285 +  if( isError ){
286 +    sprintf( painCave.errMsg,
287 +             "Error completing Shape FF in fortran\n");
288 +    painCave.isFatal = 1;
289 +    simError();
290 +  }
291 +  
292   #ifdef IS_MPI
293    sprintf( checkPointMsg,
294             "Shapes_FF atom structures successfully sent to fortran\n" );
# Line 296 | Line 297 | void SHAPES_FF::initializeAtoms( int nAtoms, Atom** th
297  
298   }
299  
300 < void SHAPES_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
300 > void Shapes_FF::cleanMe( void ){
301 >
302 > }
303 >
304 > void Shapes_FF::initializeAtoms( int nAtoms, Atom** the_atoms ){
305    
306    int i,j,k;
307 +  map<string, AtomType*>::iterator iter;
308  
309    // initialize the atoms
310    DirectionalAtom* dAtom;
311    AtomType* at;
312    DirectionalAtomType* dat;
313 <  double mySigma;
313 >  ShapeAtomType* sat;
314 >  double longCutoff;
315    double ji[3];
316    double inertialMat[3][3];
317    Mat3x3d momInt;
318    string myTypeString;
319  
320 +  bigContact = 0.0;
321 +
322    for( i=0; i<nAtoms; i++ ){
323      
324      myTypeString = the_atoms[i]->getType();
# Line 327 | Line 336 | void SHAPES_FF::initializeAtoms( int nAtoms, Atom** th
336  
337        the_atoms[i]->setMass( at->getMass() );
338        the_atoms[i]->setIdent( at->getIdent() );
339 <
340 <      if( at->isLennardJones() ) {
341 <        mySigma = at->properties->getPropertyByName("sigma");
342 <        if( bigSigma < mySigma ) bigSigma = mySigma;
339 >      
340 >      if ( at->isShape() ) {
341 >        
342 >        sat = (ShapeAtomType*)at;
343 >        longCutoff = findCutoffDistance(sat);
344 >        if (longCutoff > bigContact) bigContact = longCutoff;  
345 >        cout << bigContact << " is the cutoff value\n";
346 >        
347 >        entry_plug->useShapes = 1;
348        }
349  
350        the_atoms[i]->setHasCharge(at->isCharge());
# Line 339 | Line 353 | void SHAPES_FF::initializeAtoms( int nAtoms, Atom** th
353  
354          dat = (DirectionalAtomType*)at;
355          dAtom = (DirectionalAtom *) the_atoms[i];
356 <
356 >        
357          momInt = dat->getI();
358  
359          // zero out the moments of inertia matrix
# Line 395 | Line 409 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
409    }
410   }
411  
412 < int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAtomType* st){
412 > void Shapes_FF::parseShapeFile(string shapeFileName, ShapeAtomType* st){
413    const int MAXLEN = 1024;
414    char inLine[MAXLEN];  
415    char *token;
416    char *delim = " ,;\t\n";
417 +  int nTokens;
418    Mat3x3d momInert;
419 <  RealSphericalHarmonic tempHarmonic;
420 <  vector<RealSphericalHarmonic> functionVector;
406 <
419 >  RealSphericalHarmonic* rsh;
420 >  vector<RealSphericalHarmonic*> functionVector;
421    ifstrstream shapeFile;
422 +  string tempString;
423 +
424 +  shapeFile.open( shapeFileName.c_str() );
425 +  
426 +  if( shapeFile == NULL ){
427 +    
428 +    tempString = ffPath;
429 +    tempString += "/";
430 +    tempString += shapeFileName;
431 +    shapeFileName = tempString;
432 +        
433 +    shapeFile.open( shapeFileName.c_str() );
434 +    
435 +    if( shapeFile == NULL ){
436 +      
437 +      sprintf( painCave.errMsg,
438 +               "Error opening the shape file:\n"
439 +               "\t%s\n"
440 +               "\tHave you tried setting the FORCE_PARAM_PATH environment "
441 +               "variable?\n",
442 +               shapeFileName.c_str() );
443 +      painCave.severity = OOPSE_ERROR;
444 +      painCave.isFatal = 1;
445 +      simError();
446 +    }
447 +  }
448    
449 +  // read in the shape types.
450 +  
451    // first grab the values in the ShapeInfo section
452 <  findBegin(shapeFile, "ShapeInfo");
452 >  findBegin( shapeFile, "ShapeInfo");
453    
454    shapeFile.getline(inLine, MAXLEN);
455    while( !shapeFile.eof() ) {
# Line 416 | Line 458 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
458        // end marks section completion
459        if (isEndLine(inLine)) break;
460        
461 <      nTokens = count_tokens(inLine, delim);
461 >      nTokens = countTokens(inLine, delim);
462        if (nTokens != 0) {
463          if (nTokens < 5) {
464            sprintf( painCave.errMsg,
465 <                   "Not enough data on a ShapeInfo line in file: %s\n"
466 <                   fileName );
465 >                   "Not enough data on a ShapeInfo line in file: %s\n",
466 >                   shapeFileName.c_str() );
467            painCave.severity = OOPSE_ERROR;
468            painCave.isFatal = 1;
469            simError();
470 <
470 >        } else {
471            token = strtok(inLine, delim);
472            token = strtok(NULL, delim);
473            st->setMass(atof(token));
# Line 439 | Line 481 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
481          }
482        }
483      }
484 +    shapeFile.getline(inLine, MAXLEN);
485    }
486  
487    // now grab the contact functions
# Line 451 | Line 494 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
494      if( inLine[0] != '!' && inLine[0] != '#' ){
495        // end marks section completion
496        if (isEndLine(inLine)) break;
497 <      
455 <      nTokens = count_tokens(inLine, delim);
497 >      nTokens = countTokens(inLine, delim);
498        if (nTokens != 0) {
499          if (nTokens < 4) {
500            sprintf( painCave.errMsg,
501 <                   "Not enough data on a ContactFunctions line in file: %s\n"
502 <                   fileName );
501 >                   "Not enough data on a ContactFunctions line in file: %s\n",
502 >                   shapeFileName.c_str() );
503            painCave.severity = OOPSE_ERROR;
504            painCave.isFatal = 1;
505            simError();
506 <          
506 >        } else {
507            // read in a spherical harmonic function
508            token = strtok(inLine, delim);
509 <          tempHarmonic.setL(atoi(token));
509 >          rsh = new RealSphericalHarmonic();
510 >          rsh->setL(atoi(token));
511            token = strtok(NULL, delim);
512 <          tempHarmonic.setM(atoi(token));
512 >          rsh->setM(atoi(token));
513            token = strtok(NULL, delim);
514            if (!strcasecmp("sin",token))
515 <            tempHarmonic.makeSinFunction();
515 >            rsh->makeSinFunction();
516            else
517 <            tempHarmonic.makeCosFunction();            
517 >            rsh->makeCosFunction();            
518            token = strtok(NULL, delim);
519 <          tempHarmonic.setCoefficient(atof(token));
519 >          rsh->setCoefficient(atof(token));
520            
521 <          functionVector.push_back(tempHarmonic);
521 >          functionVector.push_back(rsh);
522          }
523        }
524      }
525 +    shapeFile.getline(inLine, MAXLEN);
526    }
527  
528    // pass contact functions to ShapeType
529 +
530    st->setContactFuncs(functionVector);
531  
532    // now grab the range functions
# Line 495 | Line 540 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
540        // end marks section completion
541        if (isEndLine(inLine)) break;
542        
543 <      nTokens = count_tokens(inLine, delim);
543 >      nTokens = countTokens(inLine, delim);
544        if (nTokens != 0) {
545          if (nTokens < 4) {
546            sprintf( painCave.errMsg,
547 <                   "Not enough data on a RangeFunctions line in file: %s\n"
548 <                   fileName );
547 >                   "Not enough data on a RangeFunctions line in file: %s\n",
548 >                   shapeFileName.c_str() );
549            painCave.severity = OOPSE_ERROR;
550            painCave.isFatal = 1;
551            simError();
552 +        } else {
553            
554            // read in a spherical harmonic function
555            token = strtok(inLine, delim);
556 <          tempHarmonic.setL(atoi(token));
556 >
557 >          rsh = new RealSphericalHarmonic();
558 >          rsh->setL(atoi(token));
559            token = strtok(NULL, delim);
560 <          tempHarmonic.setM(atoi(token));
560 >          rsh->setM(atoi(token));
561            token = strtok(NULL, delim);
562            if (!strcasecmp("sin",token))
563 <            tempHarmonic.makeSinFunction();
563 >            rsh->makeSinFunction();
564            else
565 <            tempHarmonic.makeCosFunction();            
565 >            rsh->makeCosFunction();            
566            token = strtok(NULL, delim);
567 <          tempHarmonic.setCoefficient(atof(token));
567 >          rsh->setCoefficient(atof(token));
568            
569 <          functionVector.push_back(tempHarmonic);
569 >          functionVector.push_back(rsh);
570          }
571        }
572      }
573 +    shapeFile.getline(inLine, MAXLEN);
574    }
575  
576    // pass range functions to ShapeType
# Line 538 | Line 587 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
587        // end marks section completion
588        if (isEndLine(inLine)) break;
589        
590 <      nTokens = count_tokens(inLine, delim);
590 >      nTokens = countTokens(inLine, delim);
591        if (nTokens != 0) {
592          if (nTokens < 4) {
593            sprintf( painCave.errMsg,
594 <                   "Not enough data on a StrengthFunctions line in file: %s\n"
595 <                   fileName );
594 >                   "Not enough data on a StrengthFunctions line in file: %s\n",
595 >                   shapeFileName.c_str() );
596            painCave.severity = OOPSE_ERROR;
597            painCave.isFatal = 1;
598            simError();
599 +        } else {
600            
601            // read in a spherical harmonic function
602            token = strtok(inLine, delim);
603 <          tempHarmonic.setL(atoi(token));
603 >          rsh = new RealSphericalHarmonic();
604 >          rsh->setL(atoi(token));
605            token = strtok(NULL, delim);
606 <          tempHarmonic.setM(atoi(token));
606 >          rsh->setM(atoi(token));
607            token = strtok(NULL, delim);
608            if (!strcasecmp("sin",token))
609 <            tempHarmonic.makeSinFunction();
609 >            rsh->makeSinFunction();
610            else
611 <            tempHarmonic.makeCosFunction();            
611 >            rsh->makeCosFunction();            
612            token = strtok(NULL, delim);
613 <          tempHarmonic.setCoefficient(atof(token));
613 >          rsh->setCoefficient(atof(token));
614            
615 <          functionVector.push_back(tempHarmonic);
615 >          functionVector.push_back(rsh);
616          }
617        }
618      }
619 +    shapeFile.getline(inLine, MAXLEN);
620    }
621  
622    // pass strength functions to ShapeType
# Line 572 | Line 624 | int Shapes_FF::parseShapeFile(char *shapeFile, ShapeAt
624  
625    // we're done reading from this file
626    shapeFile.close();
575 –  return 0;
627   }
628  
629 + double Shapes_FF::findLargestContactDistance(ShapeAtomType* st) {
630 +  int i, j,  nSteps;
631 +  double theta, thetaStep, thetaMin, costheta;
632 +  double phi, phiStep;
633 +  double sigma, bs;
634 +  
635 +  nSteps = 16;
636 +
637 +  thetaStep = M_PI / nSteps;
638 +  thetaMin = thetaStep / 2.0;
639 +  phiStep = thetaStep * 2.0;
640 +  bs = 0.0;
641 +  
642 +  for (i = 0; i < nSteps; i++) {
643 +    
644 +    theta = thetaMin + i * thetaStep;
645 +    costheta = cos(theta);
646 +
647 +    for (j = 0; j < nSteps; j++) {
648 +
649 +      phi = j*phiStep;
650 +
651 +      sigma = st->getContactValueAt(costheta, phi);
652 +      
653 +      if (sigma > bs) bs = sigma;
654 +    }
655 +  }
656 +
657 +  return bs;  
658 + }
659 +
660 +
661 + double Shapes_FF::findCutoffDistance(ShapeAtomType* st) {
662 +  int i, j,  nSteps;
663 +  double theta, thetaStep, thetaMin, costheta;
664 +  double phi, phiStep;
665 +  double sigma, range;
666 +  double bigCut, tempCut;
667 +
668 +  nSteps = 16;
669 +
670 +  thetaStep = M_PI / nSteps;
671 +  thetaMin = thetaStep / 2.0;
672 +  phiStep = thetaStep * 2.0;
673 +  bigCut = 0.0;
674 +  
675 +  for (i = 0; i < nSteps; i++) {
676 +    
677 +    theta = thetaMin + i * thetaStep;
678 +    costheta = cos(theta);
679 +
680 +    for (j = 0; j < nSteps; j++) {
681 +
682 +      phi = j*phiStep;
683 +
684 +      sigma = st->getContactValueAt(costheta, phi);
685 +      range = st->getRangeValueAt(costheta, phi);
686 +
687 +       // cutoff for a shape is taken to be (2.5*rangeVal + contactVal)
688 +      tempCut = 2.5*range + sigma;
689 +
690 +      if (tempCut > bigCut) bigCut = tempCut;
691 +    }
692 +  }
693 +
694 +  return bigCut;  
695 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines