~ubuntu-branches/ubuntu/vivid/gyoto/vivid

« back to all changes in this revision

Viewing changes to lib/PageThorneDisk.C

• browse files at revision 10
• compare with another revision
• download diff
• download tarball
• view history from revision 10

Show diffs side-by-side

added

removed

lib/PageThorneDisk.C

43

43

 

44

44

PageThorneDisk::PageThorneDisk() :

45

45

  ThinDisk("PageThorneDisk"), aa_(0.), aa2_(0.),

46

 

  x0_(0.), x1_(0.), x2_(0.), x3_(0.), rednoise_(0),

47

 

  uniflux_(0)

 

46

  x0_(0.), x1_(0.), x2_(0.), x3_(0.), blackbody_(0), mdot_(0),

 

47

  uniflux_(0), spectrumBB_(NULL)

48

48

{

49

49

  if (debug()) cerr << "DEBUG: PageThorneDisk Construction" << endl;

 

50

  spectrumBB_ = new Spectrum::BlackBody(); 

50

51

}

51

52

 

52

53

PageThorneDisk::PageThorneDisk(const PageThorneDisk& o) :

53

54

  ThinDisk(o), aa_(o.aa_), aa2_(o.aa2_),

54

55

  x0_(o.x0_), x1_(o.x1_), x2_(o.x2_), x3_(o.x3_),

55

 

  rednoise_(o.rednoise_), uniflux_(o.uniflux_)

 

56

  blackbody_(o.blackbody_), mdot_(0), uniflux_(o.uniflux_), 

 

57

  spectrumBB_(NULL)

56

58

{

 

59

  if (o.spectrumBB_()) spectrumBB_=o.spectrumBB_->clone();

57

60

  if (o.gg_()) gg_=o.gg_->clone();

58

61

  Generic::gg_=gg_;

59

62

  gg_->hook(this);

68

71

 

69

72

void PageThorneDisk::updateSpin() {

70

73

  if (!gg_) return;

71

 

  switch (gg_->getCoordKind()) {

 

74

  switch (gg_->coordKind()) {

72

75

  case GYOTO_COORDKIND_SPHERICAL:

73

 

    aa_ = static_cast<SmartPointer<Metric::KerrBL> >(gg_) -> getSpin();

 

76

    aa_ = static_cast<SmartPointer<Metric::KerrBL> >(gg_) -> spin();

74

77

    break;

75

78

  case GYOTO_COORDKIND_CARTESIAN:

76

 

    aa_ = static_cast<SmartPointer<Metric::KerrKS> >(gg_) -> getSpin();

 

79

    aa_ = static_cast<SmartPointer<Metric::KerrKS> >(gg_) -> spin();

77

80

    break;

78

81

  default:

79

82

    throwError("PageThorneDisk::getSpin(): unknown COORDKIND");

88

91

  x2_ = 2.*cos(acosaao3 + M_PI/3.); 

89

92

  x3_ = -2.*cos(acosaao3);

90

93

 

91

 

  rin_=(3.+z2-sqrt((3.-z1)*(3.+z1+2.*z2)));

 

94

  if (rin_==0.) rin_=(3.+z2-sqrt((3.-z1)*(3.+z1+2.*z2)));

92

95

}

93

96

 

94

 

void PageThorneDisk::setMetric(SmartPointer<Metric::Generic> gg) {

 

97

void PageThorneDisk::metric(SmartPointer<Metric::Generic> gg) {

95

98

  if (gg_) gg_->unhook(this);

96

 

  string kind = gg->getKind();

97

 

  if (kind != "KerrBL" && kind != "KerrKS" && kind != "ChernSimons")

 

99

  string kin = gg->kind();

 

100

  if (kin != "KerrBL" && kin != "KerrKS" && kin != "ChernSimons")

98

101

    throwError

99

 

      ("PageThorneDisk::setMetric(): metric must be KerrBL or KerrKS");

100

 

  ThinDisk::setMetric(gg);

 

102

      ("PageThorneDisk::metric(): metric must be KerrBL or KerrKS");

 

103

  ThinDisk::metric(gg);

101

104

  updateSpin();

102

105

  gg->hook(this);

103

106

}

105

108

double PageThorneDisk::emission(double nu_em, double dsem,

106

109

                                    double *,

107

110

                                    double coord_obj[8]) const{

108

 

  throwError("not implemented");

 

111

  if (!blackbody_) {

 

112

    throwError("In PageThorneDisk::emission: "

 

113

               "blackbody is necessary to compute emission, "

 

114

               "else, use bolometricEmission");

 

115

  }

 

116

 

 

117

  double Ibolo=bolometricEmission(nu_em,dsem,coord_obj);

 

118

  /*

 

119

    From Ibolo, find T, then Bnu(T)

 

120

   */

 

121

  double mass=gg_->mass()*1e3; // in cgs

 

122

  double c6=GYOTO_C_CGS*GYOTO_C_CGS*GYOTO_C_CGS

 

123

    *GYOTO_C_CGS*GYOTO_C_CGS*GYOTO_C_CGS;

 

124

  double g2m2=GYOTO_G_CGS*GYOTO_G_CGS*mass*mass;

 

125

  Ibolo*=mdot_*c6/g2m2; // Ibolo in cgs

 

126

  //F = sigma * T^4 (and F=pi*I)

 

127

  double TT=pow(Ibolo*M_PI/GYOTO_STEFANBOLTZMANN_CGS,0.25);

 

128

  spectrumBB_->temperature(TT);

 

129

  double Iem=(*spectrumBB_)(nu_em);

 

130

  //cout << "r T nu Iem = " << coord_obj[1] << " " << TT << " " << nu_em << " " << Iem << endl;

 

131

  if (Iem < 0.) throwError("In PageThorneDisk::emission"

 

132

                           " blackbody emission is negative!");

 

133

  return Iem;

109

134

}

110

135

 

111

136

Quantity_t PageThorneDisk::getDefaultQuantities() {

112

137

  return GYOTO_QUANTITY_USER4;

113

138

}

114

139

 

115

 

double PageThorneDisk::bolometricEmission(double dsem,

 

140

double PageThorneDisk::bolometricEmission(double /* nuem */, double dsem,

116

141

                                    double coord_obj[8]) const{

117

142

  //See Page & Thorne 74 Eqs. 11b, 14, 15. This is F(r).

118

143

  // Important remark: this emision function gives I(r),

119

144

  // not I_nu(r). And I(r)/nu^4 is conserved.

120

145

  //  cout << "r hit= " << coord_obj[1] << endl;

121

146

  if (uniflux_) return 1;

122

 

 

123

147

  double xx;

124

 

  switch (gg_->getCoordKind()) {

 

148

  switch (gg_->coordKind()) {

125

149

  case GYOTO_COORDKIND_SPHERICAL:

126

150

    xx=sqrt(coord_obj[1]);

127

151

    break;

163

187

  if (flag_radtransf_) Iem *= dsem;

164

188

  GYOTO_DEBUG_EXPR(Iem);

165

189

  

166

 

  if (rednoise_) {

167

 

    double rr=coord_obj[1], 

168

 

      r32 = pow(rr,1.5),

169

 

      tt=coord_obj[0], 

170

 

      cst_aa=0.1*sqrt(6), // amplitude cst

171

 

      cst_tt=1.;  // period cst

172

 

    srand (time(NULL));

173

 

    double random = double(rand() % 100)/100. ;

174

 

    Iem*=1.+cst_aa*random*pow(rr,-0.5)*sin(cst_tt*1./r32*tt);

175

 

    if (Iem < 0.) throwError("In PageThorneDisk::bolometricEmission"

176

 

                             " rednoised emission is negative!");

177

 

  }

178

 

  

179

190

  return Iem;

180

191

  

181

192

}

194

205

      (see below) ; this freqObs is used to transform the null

195

206

      worldline parameter dlambda (see below)

196

207

  */

197

 

  double freqObs=ph->getFreqObs(); // this is a useless quantity, always 1

 

208

  double freqObs=ph->freqObs(); // this is a useless quantity, always 1

 

209

  SmartPointer<Spectrometer::Generic> spr = ph -> spectrometer();

 

210

  size_t nbnuobs = spr() ? spr -> nSamples() : 0 ;

 

211

  double const * const nuobs = nbnuobs ? spr -> getMidpoints() : NULL;

198

212

  double dlambda = dt/coord_ph_hit[4]; //dlambda = dt/tdot

199

213

  double ggredm1 = -gg_->ScalarProd(coord_ph_hit,coord_obj_hit+4,

200

214

                                    coord_ph_hit+4);// / 1.; 

232

246

                << ") = " << dlambda << ", dsem=" << dsem << endl;

233

247

#endif

234

248

    if (data->intensity) throwError("unimplemented");

235

 

    else if (data->user4) {

236

 

      inc = (bolometricEmission(dsem, coord_obj_hit))

 

249

    if (data->user4) {

 

250

      inc = (bolometricEmission(freqObs*ggredm1, dsem, coord_obj_hit))

237

251

        * (ph -> getTransmission(size_t(-1)))

238

252

        * ggred*ggred*ggred*ggred; // I/nu^4 invariant

239

253

      *data->user4 += inc;

243

257

 

244

258

    }

245

259

    if (data->binspectrum) throwError("unimplemented");

246

 

    if (data->spectrum)  throwError("unimplemented");

 

260

    if (data->spectrum)  {

 

261

      if (!blackbody_) {

 

262

        throwError("In PageThorneDisk::process: "

 

263

                   "blackbody is necessary to compute spectrum");

 

264

      }

 

265

      for (size_t ii=0; ii<nbnuobs; ++ii) {

 

266

        double nuem=nuobs[ii]*ggredm1;

 

267

        double * tmp;

 

268

        inc = (emission(nuem, dsem, tmp, coord_obj_hit))

 

269

          * (ph -> getTransmission(size_t(-1)))

 

270

          * ggred*ggred*ggred; // Inu/nu^3 invariant

 

271

        data->spectrum[ii*data->offset] += inc;

 

272

        //cout << "in spec stored= " << ggred << " " << inc << endl;

 

273

      }

 

274

    }

247

275

    /* update photon's transmission */

248

276

    ph -> transmit(size_t(-1),

249

277

                   transmission(freqObs*ggredm1, dsem,coord_ph_hit));

255

283

}

256

284

 

257

285

void PageThorneDisk::tell(Hook::Teller* msg) {

258

 

  updateSpin();

 

286

  if (msg==gg_) updateSpin();

259

287

}

260

288

 

261

289

int PageThorneDisk::setParameter(std::string name,

262

290

                                 std::string content,

263

291

                                 std::string unit) {

264

292

  char* tc = const_cast<char*>(content.c_str());

265

 

  if (name=="RedNoise") rednoise_=1;

 

293

  if (name=="BlackbodyMdot") {

 

294

    blackbody_=1;

 

295

    mdot_=atof(tc);

 

296

  }

266

297

  if (name=="UniFlux") uniflux_=1;

267

298

  else return ThinDisk::setParameter(name, content, unit);

268

299

  return 0;

270

301

 

271

302

#ifdef GYOTO_USE_XERCES

272

303

void PageThorneDisk::fillElement(FactoryMessenger *fmp) const {

273

 

  fmp->setMetric(gg_);

 

304

  fmp->metric(gg_);

274

305

  ThinDisk::fillElement(fmp);

275

306

}

276

307

#endif

• Older »

Loggerhead is a web-based interface for Breezy
Version: 2.0.1