source: sasview/realSpaceModeling/geoshapespy/libgeoshapespy/single_helix.cc @ cd3d15b

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since cd3d15b was f2d6445, checked in by Mathieu Doucet <doucetm@…>, 17 years ago

moving realSpaceModeling to trunk

  • Property mode set to 100644
File size: 2.5 KB
Line 
1/** \file SingleHelix.cc */
2#include <cmath>
3#include <cassert>
4#include "single_helix.h"
5#include "myutil.h"
6
7SingleHelix::SingleHelix()
8{
9  hr_ = 0;
10  tr_ = 0;
11  pitch_ = 0;
12  turns_ = 0;
13}
14
15SingleHelix::SingleHelix(double helix_radius,double tube_radius,double pitch,double turns)
16{
17  hr_ = helix_radius;
18  tr_ = tube_radius;
19  pitch_ = pitch;
20  turns_ =  turns;
21}
22
23SingleHelix::~SingleHelix()
24{
25}
26
27void SingleHelix::SetHelixRadius(double hr)
28{
29  hr_ = hr;
30}
31
32void SingleHelix::SetTubeRadius(double tr)
33{
34  tr_ = tr;
35}
36
37void SingleHelix::SetPitch(double pitch)
38{
39  pitch_ = pitch;
40}
41
42void SingleHelix::SetTurns(double turns)
43{
44  turns_ = turns;
45}
46
47double SingleHelix::GetHelixRadius()
48{
49  return hr_;
50}
51
52double SingleHelix::GetTubeRadius()
53{
54  return tr_;
55}
56
57double SingleHelix::GetPitch()
58{
59  return pitch_;
60}
61
62double SingleHelix::GetTurns()
63{
64  return turns_;
65}
66
67double SingleHelix::GetMaxRadius()
68{
69  // put helix into a cylinder
70  double r_ = hr_ + tr_;
71  double l_ = pitch_*turns_ + 2*tr_;
72  double maxr = sqrt(4*r_*r_ + l_*l_)/2;
73  return maxr;
74}
75
76ShapeType SingleHelix::GetShapeType() const
77{
78  return SINGLEHELIX;
79}
80
81double SingleHelix::GetVolume()
82{
83  double V = pi*square(tr_)*sqrt(square(2*pi*hr_)+square(turns_*pitch_));
84  return V;
85}
86
87void SingleHelix::GetFormFactor(IQ * iq)
88{
89  /** number of I for output, equal to the number of rows of array IQ*/
90  /** to be finished */
91}
92
93Point3D SingleHelix::GetAPoint(double sld)
94{
95  int max_try = 100;
96  int i = 0;
97  double point1 = 0, point2 = 0, point3 = 0;
98 
99  do{
100    i++;
101    if (i > max_try)
102        break;
103    zr_ = tr_*sqrt(1+square((pitch_/(2*pi*hr_))));
104    point1 = (ran1()-0.5)*2*tr_;
105    point2 = (ran1()-0.5)*2*zr_;
106    point3 = (ran1()-0.5)*4*pi*turns_;
107  } while (((square(point1/hr_)+square(point2/zr_))>1) || (point2+pitch_*point3/(2*pi)<0));
108
109  double x = (point1 + hr_)*cos(point3);
110  double y = (point1 + hr_)*sin(point3);
111  //"-pitch_*turns_/2" is just to corretly center the helix
112  double z = point2 + pitch_*point3/(2*pi)-pitch_*turns_/2;
113
114  Point3D apoint(x,y,z,sld);
115  return apoint;
116
117 
118  std::cerr << "Max try "
119            << max_try
120            << " is reached while generating a point in cylinder" << std::endl;
121 
122  return Point3D(0,0,0);
123}
124
125bool SingleHelix::IsInside(const Point3D& point) const
126{
127  double x = point.getX()-center_[0];
128  double y = point.getY()-center_[1];
129  double z = point.getZ()-center_[2];
130
131  double p3 = atan(y/x);
132  double p2 = z-(pitch_*p3)/(2*pi);
133  double p1 = x/cos(p3)-hr_;
134
135  if ((square(p1/tr_)+square(p2/zr_))>1 || p2+pitch_*p3/(2*pi))
136    return false;
137  else
138    return true;
139}
140
Note: See TracBrowser for help on using the repository browser.