1 // qlobular.cpp
2 //
3 // Copyright (C) 2008, Celestia Development Team
4 // Initial code by Dr. Fridger Schrempp <fridger.schrempp@desy.de>
5 //
6 // Simulation of globular clusters, theoretical framework by
7 // Ivan King, Astron. J. 67 (1962) 471; ibid. 71 (1966) 64
8 //
9 // This program is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU General Public License
11 // as published by the Free Software Foundation; either version 2
12 // of the License, or (at your option) any later version.
13
14 #include <fstream>
15 #include <algorithm>
16 #include <cstdio>
17 #include <cassert>
18 #include "celestia.h"
19 #include <celmath/mathlib.h>
20 #include <celmath/perlin.h>
21 #include <celmath/intersect.h>
22 #include "astro.h"
23 #include "globular.h"
24 #include <celutil/util.h>
25 #include <celutil/debug.h>
26 #include "gl.h"
27 #include "vecgl.h"
28 #include "render.h"
29 #include "texture.h"
30 #include <math.h>
31 using namespace std;
32
33 static int cntrTexWidth = 512, cntrTexHeight = 512;
34 static int starTexWidth = 128, starTexHeight = 128;
35 static Color colorTable[256];
36 static const unsigned int GLOBULAR_POINTS = 8192;
37 static const float LumiShape = 3.0f, Lumi0 = exp(-LumiShape);
38
39 // Reference values ( = data base averages) of core radius, King concentration
40 // and mu25 isophote radius:
41
42 static const float R_c_ref = 0.83f, C_ref = 2.1f, R_mu25 = 40.32f;
43
44 // min/max c-values of globular cluster data
45
46 static const float MinC = 0.50f, MaxC = 2.58f, BinWidth = (MaxC - MinC) / 8.0f + 0.02f;
47
48 // P1 determines the zoom level, where individual cluster stars start to appear.
49 // The smaller P2 (< 1), the faster stars show up when resolution increases.
50
51 static const float P1 = 65.0f, P2 = 0.75f;
52
53 static const float RRatio_min = pow(10.0f, 1.7f);
54 static float CBin, RRatio, XI, DiskSizeInPixels, Rr = 1.0f, Gg = 1.0f, Bb = 1.0f;
55
56 static GlobularForm** globularForms = NULL;
57 static Texture* globularTex = NULL;
58 static Texture* centerTex[8] = {NULL};
59 static void InitializeForms();
60 static GlobularForm* buildGlobularForms(float);
61 static bool formsInitialized = false;
62
decreasing(const GBlob & b1,const GBlob & b2)63 static bool decreasing (const GBlob& b1, const GBlob& b2)
64 {
65 return (b1.radius_2d > b2.radius_2d);
66 }
67
GlobularTextureEval(float u,float v,float,unsigned char * pixel)68 static void GlobularTextureEval(float u, float v, float /*w*/, unsigned char *pixel)
69 {
70 // use an exponential luminosity shape for the individual stars
71 // giving sort of a halo for the brighter (i.e.bigger) stars.
72
73 float lumi = exp(- LumiShape * sqrt(u * u + v * v)) - Lumi0;
74
75 if (lumi <= 0.0f)
76 lumi = 0.0f;
77
78 int pixVal = (int) (lumi * 255.99f);
79 #ifdef HDR_COMPRESS
80 pixel[0] = 127;
81 pixel[1] = 127;
82 pixel[2] = 127;
83 #else
84 pixel[0] = 255;
85 pixel[1] = 255;
86 pixel[2] = 255;
87 #endif
88 pixel[3] = pixVal;
89 }
90
91
relStarDensity(float eta)92 float relStarDensity(float eta)
93 {
94 /*! As alpha blending weight (relStarDensity) I take the theoretical
95 * number of globular stars in 2d projection at a distance
96 * rho = r / r_c = eta * r_t from the center (cf. King_1962's Eq.(18)),
97 * divided by the area = PI * rho * rho . This number density of stars
98 * I normalized to 1 at rho=0.
99
100 * The resulting blending weight increases strongly -> 1 if the
101 * 2d number density of stars rises, i.e for rho -> 0.
102 */
103 // Since the central "cloud" is due to lack of visual resolution,
104 // rather than cluster morphology, we limit it's size by
105 // taking max(C_ref, CBin). Smaller c gives a shallower distribution!
106
107 float rRatio = max(RRatio_min, RRatio);
108 float Xi = 1.0f / sqrt(1.0f + rRatio * rRatio);
109 float XI2 = Xi * Xi;
110 float rho2 = 1.0001f + eta * eta * rRatio * rRatio; //add 1e-4 as regulator near rho=0
111
112 return ((log(rho2) + 4.0f * (1.0f - sqrt(rho2)) * Xi) / (rho2 - 1.0f) + XI2) / (1.0f - 2.0f * Xi + XI2);
113 }
114
CenterCloudTexEval(float u,float v,float,unsigned char * pixel)115 static void CenterCloudTexEval(float u, float v, float /*w*/, unsigned char *pixel)
116 {
117 /*! For reasons of speed, calculate central "cloud" texture only for
118 * 8 bins of King_1962 concentration, c = CBin, XI(CBin), RRatio(CBin).
119 */
120
121 // Skyplane projected King_1962 profile at center (rho = eta = 0):
122 float c2d = 1.0f - XI;
123
124 float eta = sqrt(u * u + v * v); // u,v = (-1..1)
125
126 // eta^2 = u * u + v * v = 1 is the biggest circle fitting into the quadratic
127 // procedural texture. Hence clipping
128
129 if (eta >= 1.0f)
130 eta = 1.0f;
131
132 // eta = 1 corresponds to tidalRadius:
133
134 float rho = eta * RRatio;
135 float rho2 = 1.0f + rho * rho;
136
137 // Skyplane projected King_1962 profile (Eq.(14)), vanishes for eta = 1:
138 // i.e. absolutely no globular stars for r > tidalRadius:
139
140 float profile_2d = (1.0f / sqrt(rho2) - 1.0f)/c2d + 1.0f ;
141 profile_2d = profile_2d * profile_2d;
142
143 #ifdef HDR_COMPRESS
144 pixel[0] = 127;
145 pixel[1] = 127;
146 pixel[2] = 127;
147 #else
148 pixel[0] = 255;
149 pixel[1] = 255;
150 pixel[2] = 255;
151 #endif
152 pixel[3] = (int) (relStarDensity(eta) * profile_2d * 255.99f);
153 }
154
Globular()155 Globular::Globular() :
156 detail (1.0f),
157 customTmpName (NULL),
158 form (NULL),
159 r_c (R_c_ref),
160 c (C_ref),
161 tidalRadius(0.0f)
162 {
163 recomputeTidalRadius();
164 }
165
cSlot(float conc) const166 unsigned int Globular::cSlot(float conc) const
167 {
168 // map the physical range of c, minC <= c <= maxC,
169 // to 8 integers (bin numbers), 0 < cSlot <= 7:
170
171 if (conc <= MinC)
172 conc = MinC;
173 if (conc >= MaxC)
174 conc = MaxC;
175
176 return (unsigned int) floor((conc - MinC) / BinWidth);
177 }
178
179
getType() const180 const char* Globular::getType() const
181 {
182 return "Globular";
183 }
184
185
setType(const std::string &)186 void Globular::setType(const std::string& /*typeStr*/)
187 {
188 }
189
getDetail() const190 float Globular::getDetail() const
191 {
192 return detail;
193 }
194
195
setDetail(float d)196 void Globular::setDetail(float d)
197 {
198 detail = d;
199 }
200
getCustomTmpName() const201 string Globular::getCustomTmpName() const
202 {
203 if (customTmpName == NULL)
204 return "";
205 else
206 return *customTmpName;
207 }
208
setCustomTmpName(const string & tmpNameStr)209 void Globular::setCustomTmpName(const string& tmpNameStr)
210 {
211 if (customTmpName == NULL)
212 customTmpName = new string(tmpNameStr);
213 else
214 *customTmpName = tmpNameStr;
215 }
216
getCoreRadius() const217 float Globular::getCoreRadius() const
218 {
219 return r_c;
220 }
221
setCoreRadius(const float coreRadius)222 void Globular::setCoreRadius(const float coreRadius)
223 {
224 r_c = coreRadius;
225 recomputeTidalRadius();
226 }
227
getHalfMassRadius() const228 float Globular::getHalfMassRadius() const
229 {
230 // Aproximation to the half-mass radius r_h [ly]
231 // (~ 20% accuracy)
232
233 return std::tan(degToRad(r_c / 60.0f)) * (float) getPosition().distanceFromOrigin() * pow(10.0f, 0.6f * c - 0.4f);
234 }
235
getConcentration() const236 float Globular::getConcentration() const
237 {
238 return c;
239 }
setConcentration(const float conc)240 void Globular::setConcentration(const float conc)
241 {
242 c = conc;
243 if (!formsInitialized)
244 InitializeForms();
245
246 // For saving time, account for the c dependence via 8 bins only,
247
248 form = globularForms[cSlot(conc)];
249 recomputeTidalRadius();
250 }
251
252
getDescription(char * buf,size_t bufLength) const253 size_t Globular::getDescription(char* buf, size_t bufLength) const
254 {
255 return snprintf(buf, bufLength, _("Globular (core radius: %4.2f', King concentration: %4.2f)"), r_c, c);
256 }
257
258
getForm() const259 GlobularForm* Globular::getForm() const
260 {
261 return form;
262 }
263
getObjTypeName() const264 const char* Globular::getObjTypeName() const
265 {
266 return "globular";
267 }
268
269
270 static const float RADIUS_CORRECTION = 0.025f;
pick(const Ray3d & ray,double & distanceToPicker,double & cosAngleToBoundCenter) const271 bool Globular::pick(const Ray3d& ray,
272 double& distanceToPicker,
273 double& cosAngleToBoundCenter) const
274 {
275 if (!isVisible())
276 return false;
277 /*
278 * The selection ellipsoid should be slightly larger to compensate for the fact
279 * that blobs are considered points when globulars are built, but have size
280 * when they are drawn.
281 */
282 Vec3d ellipsoidAxes(getRadius() * (form->scale.x + RADIUS_CORRECTION),
283 getRadius() * (form->scale.y + RADIUS_CORRECTION),
284 getRadius() * (form->scale.z + RADIUS_CORRECTION));
285
286 Quatf qf= getOrientation();
287 Quatd qd(qf.w, qf.x, qf.y, qf.z);
288
289 return testIntersection(Ray3d(Point3d() + (ray.origin - getPosition()), ray.direction) * conjugate(qd).toMatrix3(),
290 Ellipsoidd(ellipsoidAxes),
291 distanceToPicker,
292 cosAngleToBoundCenter);
293 }
load(AssociativeArray * params,const string & resPath)294 bool Globular::load(AssociativeArray* params, const string& resPath)
295 {
296 // Load the basic DSO parameters first
297
298 bool ok = DeepSkyObject::load(params, resPath);
299 if (!ok)
300 return false;
301
302 if(params->getNumber("Detail", detail))
303 setDetail((float) detail);
304
305 string customTmpName;
306 if(params->getString("CustomTemplate", customTmpName ))
307 setCustomTmpName(customTmpName);
308
309 if(params->getNumber("CoreRadius", r_c))
310 setCoreRadius(r_c);
311
312 if(params->getNumber("KingConcentration", c))
313 setConcentration(c);
314
315 return true;
316 }
317
318
render(const GLContext & context,const Vec3f & offset,const Quatf & viewerOrientation,float brightness,float pixelSize)319 void Globular::render(const GLContext& context,
320 const Vec3f& offset,
321 const Quatf& viewerOrientation,
322 float brightness,
323 float pixelSize)
324 {
325 renderGlobularPointSprites(context, offset, viewerOrientation, brightness, pixelSize);
326 }
327
328
renderGlobularPointSprites(const GLContext &,const Vec3f & offset,const Quatf & viewerOrientation,float brightness,float pixelSize)329 void Globular::renderGlobularPointSprites(const GLContext&,
330 const Vec3f& offset,
331 const Quatf& viewerOrientation,
332 float brightness,
333 float pixelSize)
334 {
335 if (form == NULL)
336 return;
337
338 float distanceToDSO = offset.length() - getRadius();
339 if (distanceToDSO < 0)
340 distanceToDSO = 0;
341
342 float minimumFeatureSize = 0.5f * pixelSize * distanceToDSO;
343
344 DiskSizeInPixels = getRadius() / minimumFeatureSize;
345
346 /*
347 * Is the globular's apparent size big enough to
348 * be noticeable on screen? If it's not, break right here to
349 * avoid all the overhead of the matrix transformations and
350 * GL state changes:
351 */
352
353 if (DiskSizeInPixels < 1.0f)
354 return;
355 /*
356 * When resolution (zoom) varies, the blended texture opacity is controlled by the
357 * factor 'pixelWeight'. At low resolution, the latter starts at 1, but tends to 0,
358 * if the resolution increases sufficiently (DiskSizeInPixels >= P1 pixels)!
359 * The smaller P2 (<1), the faster pixelWeight -> 0, for DiskSizeInPixels >= P1.
360 */
361
362 float pixelWeight = (DiskSizeInPixels >= P1)? 1.0f/(P2 + (1.0f - P2) * DiskSizeInPixels / P1): 1.0f;
363
364 // Use same 8 c-bins as in globularForms below!
365
366 unsigned int ic = cSlot(c);
367 CBin = MinC + ((float) ic + 0.5f) * BinWidth; // center value of (ic+1)th c-bin
368
369 RRatio = pow(10.0f, CBin);
370 XI = 1.0f / sqrt(1.0f + RRatio * RRatio);
371
372 if(centerTex[ic] == NULL)
373 {
374 centerTex[ic] = CreateProceduralTexture( cntrTexWidth, cntrTexHeight, GL_RGBA, CenterCloudTexEval);
375 }
376 assert(centerTex[ic] != NULL);
377
378 if (globularTex == NULL)
379 {
380 globularTex = CreateProceduralTexture( starTexWidth, starTexHeight, GL_RGBA,
381 GlobularTextureEval);
382 }
383 assert(globularTex != NULL);
384
385 glEnable (GL_BLEND);
386 glEnable (GL_TEXTURE_2D);
387 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
388
389 Mat3f viewMat = viewerOrientation.toMatrix3();
390 Vec3f v0 = Vec3f(-1, -1, 0) * viewMat;
391 Vec3f v1 = Vec3f( 1, -1, 0) * viewMat;
392 Vec3f v2 = Vec3f( 1, 1, 0) * viewMat;
393 Vec3f v3 = Vec3f(-1, 1, 0) * viewMat;
394
395 float tidalSize = 2 * tidalRadius;
396 Mat3f m =
397 Mat3f::scaling(form->scale) * getOrientation().toMatrix3() *
398 Mat3f::scaling(tidalSize);
399
400 vector<GBlob>* points = form->gblobs;
401 unsigned int nPoints =
402 (unsigned int) (points->size() * clamp(getDetail()));
403
404 /* Render central cloud sprite (centerTex). It fades away when
405 * distance from center or resolution increases sufficiently.
406 */
407
408 centerTex[ic]->bind();
409
410 float br = 2 * brightness;
411
412 glColor4f(Rr, Gg, Bb, min(br * pixelWeight, 1.0f));
413
414 glBegin(GL_QUADS);
415
416 glTexCoord2f(0, 0); glVertex(v0 * tidalSize);
417 glTexCoord2f(1, 0); glVertex(v1 * tidalSize);
418 glTexCoord2f(1, 1); glVertex(v2 * tidalSize);
419 glTexCoord2f(0, 1); glVertex(v3 * tidalSize);
420
421 glEnd();
422
423 /*! Next, render globular cluster via distinct "star" sprites (globularTex)
424 * for sufficiently large resolution and distance from center of globular.
425 *
426 * This RGBA texture fades away when resolution decreases (e.g. via automag!),
427 * or when distance from globular center decreases.
428 */
429
430
431 globularTex->bind();
432
433
434 int pow2 = 128; // Associate "Red Giants" with the 128 biggest star-sprites
435
436 float starSize = br * 0.5f; // Maximal size of star sprites -> "Red Giants"
437 float clipDistance = 100.0f; // observer distance [ly] from globular, where we
438 // start "morphing" the star-sprite sizes towards
439 // their physical values
440 glBegin(GL_QUADS);
441
442 for (unsigned int i = 0; i < nPoints; ++i)
443 {
444 GBlob b = (*points)[i];
445 Point3f p = b.position * m;
446 float eta_2d = b.radius_2d;
447
448 /*! Note that the [axis,angle] input in globulars.dsc transforms the
449 * 2d projected star distance r_2d in the globular frame to refer to the
450 * skyplane frame for each globular! That's what I need here.
451 *
452 * The [axis,angle] input will be needed anyway, when upgrading to
453 * account for ellipticities, with corresponding inclinations and
454 * position angles...
455 */
456
457
458 if ((i & pow2) != 0)
459 {
460 pow2 <<= 1;
461 starSize /= 1.25f;
462
463 if (starSize < minimumFeatureSize)
464 break;
465 }
466
467 float obsDistanceToStarRatio = (p + offset).distanceFromOrigin() / clipDistance;
468 float saveSize = starSize;
469
470 if (obsDistanceToStarRatio < 1.0f)
471 {
472 // "Morph" the star-sprite sizes at close observer distance such that
473 // the overdense globular core is dissolved upon closing in.
474
475 starSize = starSize * min(obsDistanceToStarRatio, 1.0f);
476 }
477
478 /* Colors of normal globular stars are given by color profile.
479 * Associate orange "Red Giant" stars with the largest sprite
480 * sizes (while pow2 = 128).
481 */
482
483 Color col = (pow2 < 256)? colorTable[255]: colorTable[b.colorIndex];
484 glColor4f(col.red(), col.green(), col.blue(),
485 min(br * (1.0f - pixelWeight * relStarDensity(eta_2d)), 1.0f));
486
487 glTexCoord2f(0, 0); glVertex(p + (v0 * starSize));
488 glTexCoord2f(1, 0); glVertex(p + (v1 * starSize));
489 glTexCoord2f(1, 1); glVertex(p + (v2 * starSize));
490 glTexCoord2f(0, 1); glVertex(p + (v3 * starSize));
491
492 starSize = saveSize;
493 }
494 glEnd();
495 }
496
getRenderMask() const497 unsigned int Globular::getRenderMask() const
498 {
499 return Renderer::ShowGlobulars;
500 }
501
getLabelMask() const502 unsigned int Globular::getLabelMask() const
503 {
504 return Renderer::GlobularLabels;
505 }
506
507
recomputeTidalRadius()508 void Globular::recomputeTidalRadius()
509 {
510 // Convert the core radius from arcminutes to light years
511 // Compute the tidal radius in light years
512
513 float coreRadiusLy = std::tan(degToRad(r_c / 60.0f)) * (float) getPosition().distanceFromOrigin();
514 tidalRadius = coreRadiusLy * std::pow(10.0f, c);
515 }
516
517
buildGlobularForms(float c)518 GlobularForm* buildGlobularForms(float c)
519 {
520 GBlob b;
521 vector<GBlob>* globularPoints = new vector<GBlob>;
522
523 float rRatio = pow(10.0f, c); // = r_t / r_c
524 float prob;
525 float cc = 1.0f + rRatio * rRatio;
526 unsigned int i = 0, k = 0;
527
528 // Value of King_1962 luminosity profile at center:
529
530 float prob0 = sqrt(cc) - 1.0f;
531
532 /*! Generate the globular star distribution randomly, according
533 * to the King_1962 surface density profile f(r), eq.(14).
534 *
535 * rho = r / r_c = eta r_t / r_c, 0 <= eta <= 1,
536 * coreRadius r_c, tidalRadius r_t, King concentration c = log10(r_t/r_c).
537 */
538
539 while (i < GLOBULAR_POINTS)
540 {
541 /*!
542 * Use a combination of the Inverse Transform method and
543 * Von Neumann's Acceptance-Rejection method for generating sprite stars
544 * with eta distributed according to the exact King luminosity profile.
545 *
546 * This algorithm leads to almost 100% efficiency for all values of
547 * parameters and variables!
548 */
549
550 float uu = Mathf::frand();
551
552 /* First step: eta distributed as inverse power distribution (~1/Z^2)
553 * that majorizes the exact King profile. Compute eta in terms of uniformly
554 * distributed variable uu! Normalization to 1 for eta -> 0.
555 */
556
557 float eta = tan(uu *atan(rRatio))/rRatio;
558
559 float rho = eta * rRatio;
560 float cH = 1.0f/(1.0f + rho * rho);
561 float Z = sqrt((1.0f + rho * rho)/cc); // scaling variable
562
563 // Express King_1962 profile in terms of the UNIVERSAL variable 0 < Z <= 1,
564
565 prob = (1.0f - 1.0f / Z) / prob0;
566 prob = prob * prob;
567
568 /* Second step: Use Acceptance-Rejection method (Von Neumann) for
569 * correcting the power distribution of eta into the exact,
570 * desired King form 'prob'!
571 */
572
573 k++;
574
575 if (Mathf::frand() < prob / cH)
576 {
577 /* Generate 3d points of globular cluster stars in polar coordinates:
578 * Distribution in eta (<=> r) according to King's profile.
579 * Uniform distribution on any spherical surface for given eta.
580 * Note: u = cos(phi) must be used as a stochastic variable to get uniformity in angle!
581 */
582 float u = Mathf::sfrand();
583 float theta = 2 * (float) PI * Mathf::frand();
584 float sthetu2 = sin(theta) * sqrt(1.0f - u * u);
585
586 // x,y,z points within -0.5..+0.5, as required for consistency:
587 b.position = 0.5f * Point3f(eta * sqrt(1.0f - u * u) * cos(theta), eta * sthetu2 , eta * u);
588
589 /*
590 * Note: 2d projection in x-z plane, according to Celestia's
591 * conventions! Hence...
592 */
593 b.radius_2d = eta * sqrt(1.0f - sthetu2 * sthetu2);
594
595 /* For now, implement only a generic spectrum for normal cluster
596 * stars, modelled from Hubble photo of M80.
597 * Blue Stragglers are qualitatively accounted for...
598 * assume color index poportional to Z as function of which the King profile
599 * becomes universal!
600 */
601
602 b.colorIndex = (unsigned int) (Z * 254);
603
604 globularPoints->push_back(b);
605 i++;
606 }
607 }
608
609 // Check for efficiency of sprite-star generation => close to 100 %!
610 //cout << "c = "<< c <<" i = " << i - 1 <<" k = " << k - 1 << " Efficiency: " << 100.0f * i / (float)k<<"%" << endl;
611
612 GlobularForm* globularForm = new GlobularForm();
613 globularForm->gblobs = globularPoints;
614 globularForm->scale = Vec3f(1.0f, 1.0f, 1.0f);
615
616 return globularForm;
617 }
618
InitializeForms()619 void InitializeForms()
620 {
621
622 // Build RGB color table, using hue, saturation, value as input.
623 // Hue in degrees.
624
625 // Location of hue transition and saturation peak in color index space:
626 int i0 = 36, i_satmax = 16;
627 // Width of hue transition in color index space:
628 int i_width = 3;
629
630 float sat_l = 0.08f, sat_h = 0.1f, hue_r = 27.0f, hue_b = 220.0f;
631
632 // Red Giant star color: i = 255:
633 // -------------------------------
634 // Convert hue, saturation and value to RGB
635
636 DeepSkyObject::hsv2rgb(&Rr, &Gg, &Bb, 25.0f, 0.65f, 1.0f);
637 colorTable[255] = Color(Rr, Gg, Bb);
638
639 // normal stars: i < 255, generic color profile for now, improve later
640 // --------------------------------------------------------------------
641 // Convert hue, saturation, value to RGB
642
643 for (int i = 254; i >=0; i--)
644 {
645 // simple qualitative saturation profile:
646 // i_satmax is value of i where sat = sat_h + sat_l maximal
647
648 float x = (float) i / (float) i_satmax, x2 = x ;
649 float sat = sat_l + 2 * sat_h /(x2 + 1.0f / x2);
650
651 // Fast transition from hue_r to hue_b at i = i0 within a width
652 // i_width in color index space:
653
654 float hue = hue_r + 0.5f * (hue_b - hue_r) * (std::tanh((float)(i - i0) / (float) i_width) + 1.0f);
655
656 DeepSkyObject::hsv2rgb(&Rr, &Gg, &Bb, hue, sat, 0.85f);
657 colorTable[i] = Color(Rr, Gg, Bb);
658 }
659 // Define globularForms corresponding to 8 different bins of King concentration c
660
661 globularForms = new GlobularForm*[8];
662
663 for (unsigned int ic = 0; ic <= 7; ++ic)
664 {
665 float CBin = MinC + ((float) ic + 0.5f) * BinWidth;
666 globularForms[ic] = buildGlobularForms(CBin);
667 }
668 formsInitialized = true;
669
670 }
671