# hub.darcs.net :: basilisk -> basilisk -> files

Basilisk source code

## root / src / spherical.h.page

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 /** # Spherical coordinates The default radius of the sphere is set to one. */ double Radius = 1.; /** Rather than the classical [mathematical convention](http://en.wikipedia.org/wiki/Spherical_coordinate_system), we use the geographic convention: *x* is the longitude within $[-180:180]$ degrees and *y* is the latitude within $[-90:90]$ degrees. $\Delta$ is the characteristic cell length expressed in the same units as *Radius*. */ map { x = x < -180. ? x + 360. : x > 180. ? x - 360. : x; Delta_x = Delta_y = Delta; Delta *= Radius*pi/180.; } /** On trees we need to define consistent refinement functions. The default restriction functions (averages) are already consistent (i.e. they preserve volume and length integrals). Mathematically we have $$cm = fm.y = \cos(\phi)$$ $$fm.x = 1$$ with $\phi$ the latitude in radians. This is the definition we use for the length scale factors *fm*. For the volume scale factor *cm*, more care needs to be taken to guarantee discrete volume conservation i.e. $$\sum_{children}cm_{child} = 4cm_{parent}$$ This is achieved by computing the exact area of a cell i.e. $$\Delta^2 \int^{\phi + d\phi/2}_{\phi - d \phi/2} \cos\phi d\phi = \Delta^2 (\sin(\phi + d\phi/2) - \sin(\phi - d\phi/2)) = \Delta^2 cm$$ */ #if TREE static void refine_cm_lonlat (Point point, scalar cm) { double phi = y*pi/180., dphi = Delta/(2.*Radius); double sinphi = sin(phi); fine(cm,0,0) = fine(cm,1,0) = (sinphi - sin(phi - dphi))/dphi; fine(cm,0,1) = fine(cm,1,1) = (sin(phi + dphi) - sinphi)/dphi; } static void refine_face_x_lonlat (Point point, scalar fm) { if (!is_refined(neighbor(-1))) fine(fm,0,0) = fine(fm,0,1) = 1.; if (!is_refined(neighbor(1)) && neighbor(1).neighbors) fine(fm,2,0) = fine(fm,2,1) = 1.; fine(fm,1,0) = fine(fm,1,1) = 1.; } static void refine_face_y_lonlat (Point point, scalar fm) { double phi = y*pi/180., dphi = Delta/(2.*Radius); if (!is_refined(neighbor(0,-1))) fine(fm,0,0) = fine(fm,1,0) = cos(phi - dphi); if (!is_refined(neighbor(0,1)) && neighbor(0,1).neighbors) fine(fm,0,2) = fine(fm,1,2) = cos(phi + dphi); fine(fm,0,1) = fine(fm,1,1) = cos(phi); } #endif event metric (i = 0) { /** We initialise the scale factors, taking care to first allocate the fields if they are still constant. */ if (is_constant(cm)) { scalar * l = list_copy (all); cm = new scalar; free (all); all = list_concat ({cm}, l); free (l); } scalar cmv = cm; foreach() { double phi = y*pi/180., dphi = Delta/(2.*Radius); cmv[] = (sin(phi + dphi) - sin(phi - dphi))/(2.*dphi); } if (is_constant(fm.x)) { scalar * l = list_copy (all); fm = new face vector; free (all); all = list_concat ((scalar *){fm}, l); free (l); } face vector fmv = fm; foreach_face(x) fmv.x[] = 1.; foreach_face(y) fmv.y[] = cos(y*pi/180.); /** We set our refinement/prolongation functions on trees. */ #if TREE cm.refine = cm.prolongation = refine_cm_lonlat; fm.x.prolongation = refine_face_x_lonlat; fm.y.prolongation = refine_face_y_lonlat; #endif boundary ({cm, fm}); }