ARTS  2.2.66
m_zeeman.cc
Go to the documentation of this file.
1 /* Copyright (C) 2012
2  Richard Larsson <ric.larsson@gmail.com>
3 
4  This program is free software; you can redistribute it and/or modify it
5  under the terms of the GNU General Public License as published by the
6  Free Software Foundation; either version 2, or (at your option) any
7  later version.
8 
9  This program is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU General Public License for more details.
13 
14  You should have received a copy of the GNU General Public License
15  along with this program; if not, write to the Free Software
16  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17  USA. */
18 
19 #include <cmath>
20 #include <stdexcept>
21 #include <string>
22 #include "auto_md.h"
23 #include "ppath.h"
24 #include "messages.h"
25 #include "math_funcs.h"
26 #include "absorption.h"
27 #include "abs_species_tags.h"
28 #include "physics_funcs.h"
29 #include "matpackIII.h"
30 #include "rte.h"
31 #include "rational.h"
32 
33 extern const Numeric PI;
34 extern const Numeric DEG2RAD;
35 extern const Numeric RAD2DEG;
36 extern const Numeric PLANCK_CONST;
37 extern const Numeric BOHR_MAGNETON;
38 extern const Numeric LANDE_GS;
39 
40 
55 void phase_matrix(MatrixView K, const Numeric& theta, const Numeric& eta, const Index& DM)
56 {
57  assert(K.nrows() == 4 );
58  assert(K.ncols() == 4 );
59 
60  const Numeric
61  S2T = sin(theta)*sin(theta),
62  CT = cos(theta),
63  CE2 = cos(2*eta),
64  SE2 = sin(2*eta);
65 
66 
67  switch( DM )
68  {
69  case -1: // Transitions anti-parallel to the magnetic field
70  K(0,0) = 0; K(0,1) = 0; K(0,2) = 0; K(0,3) = 0;
71  K(1,0) = 0; K(1,1) = 0; K(1,2) = 2 * CT; K(1,3) = S2T * SE2;
72  K(2,0) = 0; K(2,1) = - 2 * CT; K(2,2) = 0; K(2,3) = - S2T * CE2;
73  K(3,0) = 0; K(3,1) = - S2T * SE2; K(3,2) = S2T * CE2; K(3,3) = 0;
74  break;
75  case 1: // Transitions parallel to the magnetic field
76  K(0,0) = 0; K(0,1) = 0; K(0,2) = 0; K(0,3) = 0;
77  K(1,0) = 0; K(1,1) = 0; K(1,2) = - 2 * CT; K(1,3) = S2T * SE2;
78  K(2,0) = 0; K(2,1) = 2 * CT; K(2,2) = 0; K(2,3) = - S2T * CE2;
79  K(3,0) = 0; K(3,1) = - S2T * SE2; K(3,2) = S2T * CE2; K(3,3) = 0;
80  break;
81  case 0:// Transitions perpendicular to the magnetic field
82  K(0,0) = 0; K(0,1) = 0; K(0,2) = 0; K(0,3) = 0;
83  K(1,0) = 0; K(1,1) = 0; K(1,2) = 0; K(1,3) = - S2T * SE2;
84  K(2,0) = 0; K(2,1) = 0; K(2,2) = 0; K(2,3) = S2T * CE2;
85  K(3,0) = 0; K(3,1) = S2T * SE2; K(3,2) = - S2T * CE2; K(3,3) = 0;
86  break;
87  default: // Nil matrix since this should not be called.
88  K=0;
89  break;
90  };
91 };
92 
93 
108 void attenuation_matrix(MatrixView K, const Numeric& theta, const Numeric& eta, const Index& DM)
109 {
110  assert(K.nrows() == 4 );
111  assert(K.ncols() == 4 );
112 
113 
114  const Numeric
115  S2T = sin(theta)*sin(theta),
116  C2T = cos(theta)*cos(theta),
117  CT = cos(theta),
118  CE2 = cos(2*eta),
119  SE2 = sin(2*eta);
120 
121  switch( DM )
122  {
123  case -1: // Transitions anti-parallel to the magnetic field
124  K(0,0) = 1 + C2T; K(0,1) = S2T * CE2; K(0,2) = S2T * SE2; K(0,3) = 2 * CT;
125  K(1,0) = S2T * CE2; K(1,1) = 1 + C2T; K(1,2) = 0; K(1,3) = 0;
126  K(2,0) = S2T * SE2; K(2,1) = 0; K(2,2) = 1 + C2T; K(2,3) = 0;
127  K(3,0) = 2 * CT; K(3,1) = 0; K(3,2) = 0; K(3,3) = 1 + C2T;
128  break;
129  case 1: // Transitions parallel to the magnetic field
130  K(0,0) = 1 + C2T; K(0,1) = S2T * CE2; K(0,2) = S2T * SE2; K(0,3) = - 2 * CT;
131  K(1,0) = S2T * CE2; K(1,1) = 1 + C2T; K(1,2) = 0; K(1,3) = 0;
132  K(2,0) = S2T * SE2; K(2,1) = 0; K(2,2) = 1 + C2T; K(2,3) = 0;
133  K(3,0) = - 2 * CT; K(3,1) = 0; K(3,2) = 0; K(3,3) = 1 + C2T;
134  break;
135  case 0: // Transitions perpendicular to the magnetic field
136  K(0,0) = S2T; K(0,1) = - S2T * CE2; K(0,2) = - S2T * SE2; K(0,3) = 0;
137  K(1,0) = - S2T * CE2; K(1,1) = S2T; K(1,2) = 0; K(1,3) = 0;
138  K(2,0) = - S2T * SE2; K(2,1) = 0; K(2,2) = S2T; K(2,3) = 0;
139  K(3,0) = 0; K(3,1) = 0; K(3,2) = 0; K(3,3) = S2T;
140  break;
141  default: // Unity matrix in attenuation
142  K(0,0) = 1; K(0,1) = 0; K(0,2) = 0; K(0,3) = 0;
143  K(1,0) = 0; K(1,1) = 1; K(1,2) = 0; K(1,3) = 0;
144  K(2,0) = 0; K(2,1) = 0; K(2,2) = 1; K(2,3) = 0;
145  K(3,0) = 0; K(3,1) = 0; K(3,2) = 0; K(3,3) = 1;
146  break;
147  };
148 };
149 
150 
151 Numeric gs_caseb(const Numeric& N, const Numeric& J, const Numeric& S, const Numeric& GS) { return (GS*(J*(J+1.)+S*(S+1.)-N*(N+1.))/(J*(J+1.))/2.0); }
152 Numeric gs_casea(const Numeric& Omega, const Numeric& J, const Numeric& Sigma, const Numeric& GS) { return GS/2.0/J/(J+1)*Omega*(Omega+Sigma); }
153 
154 
169 Numeric relative_strength(const Rational& m, const Rational& j, const Index& dj, const Index& dm)
170 {
171  // Numeric conversions are necessary.
172  const Numeric J = (j-dj).toNumeric(), M = (m).toNumeric();
173 
174  // Variable to be returned.
175  Numeric ret_val;
176 
177  switch ( dj )
178  {
179  case -1:
180  switch ( dm )
181  {
182  case -1: // Transitions anti-parallel to the magnetic field
183  ret_val = (0.75)*(J+M)*(J-1+M)/(2.*J*(2.*J-1)*(2.*J+1));
184  break;
185  case 0: // Transitions perpendicular to the magnetic field
186  ret_val = (1.50)*(J*J-M*M)/(J*(2.*J-1.)*(2.*J+1.));
187  break;
188  case +1: // Transitions parallel to the magnetic field
189  ret_val = (0.75)*(J-M)*(J-1.-M)/(2.*J*(2.*J-1.)*(2.*J+1.));
190  break;
191  default:
192  throw std::runtime_error("Something is extremely wrong.");
193  break;
194  }
195  break;
196  case 0:
197  switch ( dm )
198  {
199  case -1: // Transitions anti-parallel to the magnetic field
200  ret_val = (0.75)*(J+M)*(J+1.-M)/(2.*J*(J+1.)*(2.*J+1.));
201  break;
202  case 0: // Transitions perpendicular to the magnetic field
203  ret_val = (1.50)*M*M/(J*(J+1.)*(2.*J+1.));
204  break;
205  case +1: // Transitions parallel to the magnetic field
206  ret_val = (0.75)*(J-M)*(J+1.+M)/(2.*J*(J+1.)*(2.*J+1.));
207  break;
208  default:
209  throw std::runtime_error("Something is extremely wrong.");
210  break;
211  }
212  break;
213  case +1:
214  switch ( dm )
215  {
216  case -1: // Transitions anti-parallel to the magnetic field
217  ret_val = (0.75)*(J+1.-M)*(J+2.-M)/(2.*(J+1.)*(2.*J+1.)*(2.*J+3.));
218  break;
219  case 0: // Transitions perpendicular to the magnetic field
220  ret_val = (1.5)*((J+1.)*(J+1.)-M*M)/((J+1.)*(2.*J+1.)*(2.*J+3.));
221  break;
222  case +1: // Transitions parallel to the magnetic field
223  ret_val = (0.75)*(J+1.+M)*(J+2.+M)/(2.*(J+1.)*(2.*J+1.)*(2.*J+3.));
224  break;
225  default:
226  throw std::runtime_error("Something is extremely wrong.");
227  break;
228  }
229  break;
230  default:
231  throw std::runtime_error("Something is extremely wrong.");
232  break;
233  }
234 
235  return ret_val;
236 }
237 
238 
239 Numeric (*frequency_change)(const Rational&, const Rational&, const Rational&, const Numeric&, const Index&, const Index&, const Index&, const Numeric&, const Numeric&);
240 
241 
261 Numeric frequency_change_caseb(const Rational& n, const Rational& m, const Rational& j, const Numeric& S, const Index& DJ, const Index& DM, const Index& DN, const Numeric& H_mag, const Numeric& GS)
262 {
263  const Numeric N_up = n.toNumeric();
264  const Numeric N_lo = N_up - (Numeric)DN;
265  const Numeric M_lo = m.toNumeric();
266  const Numeric M_up = M_lo + (Numeric)DM; // This will probably confuse people even though it is correct. Should fix so outer loop is over M_up.
267  const Numeric J_up = j.toNumeric();
268  const Numeric J_lo = J_up - (Numeric)DJ;
269 
270  // The special case is because g_s is undefined for J=0.
271  return (!(j == 0 && DJ == -1)) ?
272  - H_mag * (M_lo) * gs_caseb(N_lo,J_lo,S,GS) / PLANCK_CONST * BOHR_MAGNETON
273  + H_mag * (M_up) * gs_caseb(N_up,J_up,S,GS) / PLANCK_CONST * BOHR_MAGNETON :
274  H_mag * (M_lo) * gs_caseb(N_lo,J_lo,S,GS) / PLANCK_CONST * BOHR_MAGNETON ;
275 }
276 
277 
297 Numeric frequency_change_casea(const Rational& omega, const Rational& m, const Rational& j, const Numeric& Sigma, const Index& DJ, const Index& DM, const Index& Domega, const Numeric& H_mag, const Numeric& GS)
298 {
299  const Numeric Omega_up = omega.toNumeric();
300  const Numeric Omega_lo = Omega_up - (Numeric)Domega;
301  const Numeric M_lo = m.toNumeric();
302  const Numeric M_up = M_lo + (Numeric)DM; // This will probably confuse people even though it is correct. Should fix so outer loop is over M_up.
303  const Numeric J_up = j.toNumeric();
304  const Numeric J_lo = J_up - (Numeric)DJ;
305 
306  // This follows Berdyugina and Solnaki
307  return
308  H_mag * (M_up) * gs_casea(Omega_up,J_up,Sigma,GS) / PLANCK_CONST * BOHR_MAGNETON -
309  H_mag * (M_lo) * gs_casea(Omega_lo,J_lo,Sigma,GS) / PLANCK_CONST * BOHR_MAGNETON ;
310 }
311 
312 
316 void xsec_species_line_mixing_wrapper_with_zeeman( Tensor3View part_abs_mat, const ArrayOfArrayOfSpeciesTag& abs_species, const ArrayOfLineshapeSpec& abs_lineshape,
317  const ArrayOfLineRecord& lr, const SpeciesAuxData& isotopologue_ratios, const Matrix& abs_vmrs, const Vector& abs_p,
318  const Vector& abs_t, const Vector& f_grid, const Numeric& theta, const Numeric& eta, const Index& DM, const Index& this_species,
319  const ArrayOfArrayOfLineMixingRecord& line_mixing_data,
320  const ArrayOfArrayOfIndex& line_mixing_data_lut,
321  const ArrayOfIndex& temp_species_lut,
322  const Verbosity& verbosity )
323 {
324  assert( part_abs_mat.npages() == f_grid.nelem() && part_abs_mat.ncols() == 4 && part_abs_mat.nrows() == 4 );
325 
326  Matrix A(f_grid.nelem(), 1, 0.);
327  Matrix B(f_grid.nelem(), 1, 0.);
328 
329  // Handle line-mixing, if line-mixing is supposed to be handled.
330  ArrayOfArrayOfIndex temp_lut = line_mixing_data_lut;
331  for ( Index i=0; i<abs_species[this_species].nelem(); ++i )
332  if( abs_species[this_species][i].LineMixingType() != SpeciesTag::LINE_MIXING_TYPE_NONE )
333  { temp_lut[this_species] = temp_species_lut; break;}
334 
335  for ( Index i=0; i<abs_species[this_species].nelem(); ++i )
336  {
337  Matrix attenuation(f_grid.nelem(), 1, 0.), phase(f_grid.nelem(), 1, 0.);;
338 
339  xsec_species_line_mixing_wrapper( attenuation, phase, line_mixing_data, temp_lut,
340  f_grid, abs_p, abs_t, abs_vmrs, abs_species, this_species, lr,
341  abs_lineshape[i].Ind_ls(), abs_lineshape[i].Ind_lsn(), abs_lineshape[i].Cutoff(),
342  isotopologue_ratios, verbosity ); // Now in cross section
343 
344  attenuation *= abs_vmrs(this_species, 0) * number_density( abs_p[0],abs_t[0]); // Now in absorption coef.
345  phase *= abs_vmrs(this_species, 0) * number_density( abs_p[0],abs_t[0]); // Now in absorption coef.
346  phase *= 2; // phase matrix is twice as large according to sources.
347 
348  A += attenuation;
349  B += phase;
350  }
351  Matrix K_a(4,4), K_b(4,4);
352  attenuation_matrix(K_a, theta*DEG2RAD, eta*DEG2RAD, DM);
353  phase_matrix(K_b, theta*DEG2RAD, eta*DEG2RAD, DM);
354 
355  Tensor3 temp_part_abs_mat=part_abs_mat;
356  mult(part_abs_mat,A(joker,0), K_a);// Resets part_abs_mat
357  mult(temp_part_abs_mat,B(joker,0), K_b);
358 
359  part_abs_mat+=temp_part_abs_mat;
360 
361 }
362 
363 
364 /* Workspace method: Doxygen documentation will be auto-generated */
365 void propmat_clearskyAddZeeman(Tensor4& propmat_clearsky,
366  const Vector& f_grid,
367  const ArrayOfArrayOfSpeciesTag& abs_species,
368  const ArrayOfArrayOfLineRecord& abs_lines_per_species,
369  const ArrayOfLineshapeSpec& abs_lineshape,
370  const SpeciesAuxData& isotopologue_ratios,
371  const SpeciesAuxData& isotopologue_quantum,
372  const Numeric& rtp_pressure,
373  const Numeric& rtp_temperature,
374  const Vector& rtp_vmr,
375  const Vector& rtp_mag,
376  const Vector& ppath_los,
377  const Index& atmosphere_dim,
378  const ArrayOfArrayOfLineMixingRecord& line_mixing_data,
379  const ArrayOfArrayOfIndex& line_mixing_data_lut,
380  const Index& manual_zeeman_tag,
381  const Numeric& manual_zeeman_magnetic_field_strength,
382  const Numeric& manual_zeeman_theta,
383  const Numeric& manual_zeeman_eta,
384  const Verbosity& verbosity)
385 {
386  CREATE_OUT3;
387 
388  const Numeric margin = 1e-4; // This margin can perhaps be lowered? Perhaps by returning RS as Rational?
389  bool do_zeeman = false;
390 
391  // Check that correct isotopologue ratios are defined for the species
392  // we want to calculate
393  checkIsotopologueRatios(abs_species, isotopologue_ratios);
394 
395  {// Begin TEST(s)
396  if (abs_species.nelem() != abs_lines_per_species.nelem())
397  throw std::runtime_error("Dimension of *abs_species* and *abs_lines_per_species* don't match.");
398  for(Index II = 0; II<abs_species.nelem(); II++)
399  if(is_zeeman(abs_species[II])) { do_zeeman = true; break; } // If any species is Zeeman, do it.
400  if( propmat_clearsky.ncols() != 4 )
401  throw std::runtime_error("Zeeman Effect is only implemented for Stokes dimension 4.");
402  if( propmat_clearsky.nrows() != 4 )
403  throw std::runtime_error("Zeeman Effect is only implemented for Stokes dimension 4.");
404  if( propmat_clearsky.npages() != f_grid.nelem() )
405  throw std::runtime_error("Frequency dimension of *propmat_clearsky* not equal to length of *f_grid*.");
406  if( propmat_clearsky.nbooks() != abs_species.nelem() )
407  throw std::runtime_error("Species dimension of *propmat_clearsky* not equal to length of *abs_species*.");
408  if( rtp_mag.nelem() != 3 )
409  throw std::runtime_error("*rtp_mag* must have length 3.");
410  if( atmosphere_dim != 3 )
411  throw std::runtime_error("*atmosphere_dim* must be 3.");
412  if( ppath_los.nelem() != 2 )
413  throw std::runtime_error("*ppath_los* is not set correctly.");
414  }// End TEST(s)
415 
416  Vector R_path_los;
417  mirror_los(R_path_los, ppath_los, atmosphere_dim);
418 
419  // Using the standard scalar absorption functions to get physics parameters,
420  Vector abs_p, abs_t; Matrix abs_vmrs;
421  AbsInputFromRteScalars( abs_p, abs_t, abs_vmrs, // Output
422  rtp_pressure, rtp_temperature, rtp_vmr, //Input
423  verbosity); // Verbose!
424  if( do_zeeman && (( rtp_mag[0]!=0 || rtp_mag[1]!=0 || rtp_mag[2]!=0 ) || manual_zeeman_tag != 0 ))
425  {
426  //Get the magnitude of the magnetic field and store a local unit Vector for simplified angle calculations.
427  const Numeric H_mag = manual_zeeman_tag != 0?manual_zeeman_magnetic_field_strength:sqrt( rtp_mag * rtp_mag );
428 
429  Numeric theta, eta;
430  if(manual_zeeman_tag!=0)
431  { // Leaving it up to the user to manually tag the angles for simplified magnetic fields.
432  eta = manual_zeeman_eta;
433  theta = manual_zeeman_theta;
434  }
435  else
436  { // Getting angles from coordinate system
437  Vector H(3);
438  H = rtp_mag;
439  H /= H_mag;
440  // Direction vector of radiation
441  Numeric dx, dy, dz;
442  zaaa2cart(dx,dy,dz,R_path_los[0],R_path_los[1]);
443  // Radiation path direction as per Mishchenko.
444  Vector R_path(3);
445  R_path[0] = dx;
446  R_path[1] = dy;
447  R_path[2] = dz;
448  // Vertical polarization direction as per Mishchenko.
449  Vector e_v(3);
450  e_v[0] = cos(R_path_los[0] * DEG2RAD) * cos(R_path_los[1] * DEG2RAD);
451  e_v[1] = cos(R_path_los[0] * DEG2RAD) * sin(R_path_los[1] * DEG2RAD);
452  e_v[2] = -sin(R_path_los[0] * DEG2RAD);
453  // Horizontal polarization direction as per Mishchenko.
454  Vector e_h(3);
455  e_h[0] = -sin(R_path_los[1] * DEG2RAD);
456  e_h[1] = cos(R_path_los[1] * DEG2RAD);
457  e_h[2] = 0;
458  // Get the coordinate system used by Lenoir.
459  Vector tmp(3);
460  proj(tmp, R_path, H);
461  Vector R11(3);
462  R11 = H;
463  R11 -= tmp;
464  R11 *= sqrt(R11*R11);
465  Vector R22(3);
466  cross3(R22, R11, R_path);
467  // Test if the rotation is clockwise or counterclockwise.
468  const Numeric eta_test = vector_angle(R22, e_h);
469  // Find the angle between Mishchenko vertical/horizontal and Lenoir vertical/horizontal
470  eta = (eta_test > 90.0)?-vector_angle(R22, e_v):vector_angle(R22, e_v);
471 
472  // Angle between radiation propagation and magnetic field for determining how the radiation is polarized..
473  theta = vector_angle(H, R_path);
474  }
475 
476  Numeric DF, RS; // Delta Frequency and Relative Strength
477 
478  // For all species
479  for(Index II = 0; II<abs_species.nelem(); II++)
480  {
481  // Reinitialize every loop to empty the set.
482  ArrayOfLineRecord temp_abs_lines_sm, temp_abs_lines_sp, //sigma minus, sigma plus
483  temp_abs_lines_pi; // pi
484  ArrayOfIndex temp_lut_sm, temp_lut_sp, temp_lut_pi;
485 
486  // If the species isn't Zeeman, look at the next species
487  if(!is_zeeman(abs_species[II])) continue;
488 
489  // Else loop over all the lines in the species.
490  for (Index ii = 0; ii< abs_lines_per_species[II].nelem(); ii++)
491  {
492  // local LineRecord
493  LineRecord temp_LR = abs_lines_per_species[II][ii];
494  const Rational J = temp_LR.QuantumNumbers().Lower(QN_J);
495  const Numeric hund = isotopologue_quantum.getParam(temp_LR.Species(), temp_LR.Isotopologue(), 2);
496  Numeric S;
497  const Numeric GS = isotopologue_quantum.getParam(temp_LR.Species(), temp_LR.Isotopologue(), 0);
498  const Index DJ = (J - temp_LR.QuantumNumbers().Upper(QN_J)).toIndex();
499  Numeric RS_sum = 0; //Sum relative strength (which ought be close to one by the end)
500  Index number_M_sm = 0, number_M_sp=0, number_M_pi=0;
501  // Only look at lines which have no change in the main rotational number
502  // std::cout << "RSpi=[];RSsp=[];RSsm=[];DFpi=[];DFsm=[];DFsp=[];\n";
503 
504  Rational Main;
505  Index DMain;
506 
507  if( hund ==0 )//Case a
508  {
509  Main = temp_LR.QuantumNumbers().Lower(QN_Omega);
510  DMain = (Main - temp_LR.QuantumNumbers().Upper(QN_Omega)).toIndex();
511  frequency_change=frequency_change_casea;
512  S = 1-Main.toNumeric();
513  }
514  else if( hund == 1 )// Case b
515  {
516  Main = temp_LR.QuantumNumbers().Lower(QN_N);
517  DMain = (Main - temp_LR.QuantumNumbers().Upper(QN_N)).toIndex();
518  frequency_change=frequency_change_caseb;
519  S = isotopologue_quantum.getParam(temp_LR.Species(), temp_LR.Isotopologue(), 1);
520  }
521  else
522  {
523  ostringstream os;
524  os << "There are undefined Hund cases: " << abs_lines_per_species[II][ii] << "\nThe case is: "<<hund<<", allowed are (a): "<<0<<" and (b): " << 1<<"\n";
525  throw std::runtime_error(os.str());
526  }
527 
528  if (!J.isUndefined() != 0 && !Main.isUndefined() != 0 ) // This means the lines are considered erroneous.
529  {
530 
531  for ( Rational M = -J+DJ; M<=J-DJ; M++ )
532  {
533  /*
534  Note that:
535  sp := sigma plus, which means DM = 1
536  sm := sigma minus, which means DM = -1
537  pi := planar, which means DM = 0
538  */
539  if ( DJ == 1 )
540  { // Then all DM transitions possible for all M
541  DF = frequency_change(Main, M, J, S, DJ, -1, DMain, H_mag, GS);
542  RS = relative_strength(M, J, DJ, -1);
543  RS_sum += RS;
544  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
545  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
546  temp_abs_lines_sm.push_back(temp_LR);
547  number_M_sm++;
548  // std::cout << "RSsm=[RSsm " << RS << "];DFsm=[DFsm " << DF << "];\n";
549 
550  DF = frequency_change(Main, M, J, S, DJ, 0, DMain, H_mag, GS);
551  RS = relative_strength(M, J, DJ, 0);
552  RS_sum += RS;
553  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
554  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
555  temp_abs_lines_pi.push_back(temp_LR);
556  number_M_pi++;
557  // std::cout << "RSpi=[RSpi " << RS << "];DFpi=[DFpi " << DF << "];\n";
558 
559  DF = frequency_change(Main, M, J, S, DJ, +1, DMain, H_mag, GS);
560  RS = relative_strength(M, J, DJ, +1);
561  RS_sum += RS;
562  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
563  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
564  temp_abs_lines_sp.push_back(temp_LR);
565  number_M_sp++;
566  // std::cout << "RSsp=[RSsp " << RS << "];DFsp=[DFsp " << DF << "];\n";
567  }
568  else if ( DJ == 0 )
569  { // Then all DM transitions possible for all M
570  DF = frequency_change(Main, M, J, S, DJ, -1, DMain, H_mag, GS);
571  RS = relative_strength(M, J, DJ, -1);
572  RS_sum += RS;
573  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
574  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
575  temp_abs_lines_sm.push_back(temp_LR);
576  number_M_sm++;
577  // std::cout << "RSsm=[RSsm " << RS << "];DFsm=[DFsm " << DF << "];\n";
578  if( ! (M == 0) )
579  {
580  DF = frequency_change(Main, M, J, S, DJ, 0, DMain, H_mag, GS);
581  RS = relative_strength(M, J, DJ, 0);
582  RS_sum += RS;
583  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
584  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
585  temp_abs_lines_pi.push_back(temp_LR);
586  number_M_pi++;
587  // std::cout << "RSpi=[RSpi " << RS << "];DFpi=[DFpi " << DF << "];\n";
588  }
589 
590  DF = frequency_change(Main, M, J, S, DJ, +1, DMain, H_mag, GS);
591  RS = relative_strength(M, J, DJ, +1);
592  RS_sum += RS;
593  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
594  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
595  temp_abs_lines_sp.push_back(temp_LR);
596  number_M_sp++;
597  // std::cout << "RSsp=[RSsp " << RS << "];DFsp=[DFsp " << DF << "];\n";
598  }
599  else if ( DJ == -1 )
600  { // Then certain M results in blocked DM transitions
601  if ( M == -J + DJ && M!=0 )
602  { // Lower limit M only allows DM = 1
603  DF = frequency_change(Main, M, J, S, DJ, +1, DMain, H_mag, GS);
604  RS = relative_strength(M, J, DJ, +1);
605  RS_sum += RS;
606  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
607  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
608  temp_abs_lines_sp.push_back(temp_LR);
609  number_M_sp++;
610  // std::cout << "RSsp=[RSsp " << RS << "];DFsp=[DFsp " << DF << "];\n";
611 
612  }
613  else if ( M == -J + DJ + 1 && M!=0 )
614  { // Next to lower limit M can only allow DM = 1, 0
615  DF = frequency_change(Main, M, J, S, DJ, +1, DMain, H_mag, GS);
616  RS = relative_strength(M, J, DJ, +1);
617  RS_sum += RS;
618  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
619  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
620  temp_abs_lines_sp.push_back(temp_LR);
621  number_M_sp++;
622  // std::cout << "RSsp=[RSsp " << RS << "];DFsp=[DFsp " << DF << "];\n";
623 
624  DF = frequency_change(Main, M, J, S, DJ, 0, DMain, H_mag, GS);
625  RS = relative_strength(M, J, DJ, 0);
626  RS_sum += RS;
627  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
628  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
629  temp_abs_lines_pi.push_back(temp_LR);
630  number_M_pi++;
631  // std::cout << "RSpi=[RSpi " << RS << "];DFpi=[DFpi " << DF << "];\n";
632  }
633  else if ( M == J - DJ - 1 && M!=0 )
634  { // Next to upper limit M can only allow DM = 0, -1
635  DF = frequency_change(Main, M, J, S, DJ, 0, DMain, H_mag, GS);
636  RS = relative_strength(M, J, DJ, 0);
637  RS_sum += RS;
638  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
639  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
640  temp_abs_lines_pi.push_back(temp_LR);
641  number_M_pi++;
642  // std::cout << "RSpi=[RSpi " << RS << "];DFpi=[DFpi " << DF << "];\n";
643 
644  DF = frequency_change(Main, M, J, S, DJ, -1, DMain, H_mag, GS);
645  RS = relative_strength(M, J, DJ, -1);
646  RS_sum += RS;
647  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
648  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
649  temp_abs_lines_sm.push_back(temp_LR);
650  number_M_sm++;
651  // std::cout << "RSsm=[RSsm " << RS << "];DFsm=[DFsm " << DF << "];\n";
652  }
653  else if ( M == J - DJ && M!=0 )
654  { // Upper limit M only allow DM = -1
655  DF = frequency_change(Main, M, J, S, DJ, -1, DMain, H_mag, GS);
656  RS = relative_strength(M, J, DJ, -1);
657  RS_sum += RS;
658  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
659  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
660  temp_abs_lines_sm.push_back(temp_LR);
661  number_M_sm++;
662  // std::cout << "RSsm=[RSsm " << RS << "];DFsm=[DFsm " << DF << "];\n";
663  }
664  else if( (-J + DJ + 1) == (J - DJ - 1) && M == 0)
665  { // Special case for N=1, J=0, M=0. Only allows DM = 0
666  DF = frequency_change(Main, M, J, S, DJ, 0, DMain, H_mag, GS);
667  RS = relative_strength(M, J, DJ, 0);
668  RS_sum += RS;
669  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
670  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
671  temp_abs_lines_pi.push_back(temp_LR);
672  number_M_pi++;
673  // std::cout << "RSpi=[RSpi " << RS << "];DFpi=[DFpi " << DF << "];\n";
674  }
675  else
676  { // All DM transitions are possible for these M(s)
677  DF = frequency_change(Main, M, J, S, DJ, +1, DMain, H_mag, GS);
678  RS = relative_strength(M, J, DJ, +1);
679  RS_sum += RS;
680  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
681  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
682  temp_abs_lines_sp.push_back(temp_LR);
683  number_M_sp++;
684  // std::cout << "RSsp=[RSsp " << RS << "];DFsp=[DFsp " << DF << "];\n";
685 
686  DF = frequency_change(Main, M, J, S, DJ, 0, DMain, H_mag, GS);
687  RS = relative_strength(M, J, DJ, 0);
688  RS_sum += RS;
689  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
690  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
691  temp_abs_lines_pi.push_back(temp_LR);
692  number_M_pi++;
693  // std::cout << "RSpi=[RSpi " << RS << "];DFpi=[DFpi " << DF << "];\n";
694 
695  DF = frequency_change(Main, M, J, S, DJ, -1, DMain, H_mag, GS);
696  RS = relative_strength(M, J, DJ, -1);
697  RS_sum += RS;
698  temp_LR.setF( abs_lines_per_species[II][ii].F() + DF );
699  temp_LR.setI0( abs_lines_per_species[II][ii].I0() * RS );
700  temp_abs_lines_sm.push_back(temp_LR);
701  number_M_sm++;
702  // std::cout << "RSsm=[RSsm " << RS << "];DFsm=[DFsm " << DF << "];\n";
703  }
704  }
705  else
706  { // The tests above failed and catastrophe follows
707  ostringstream os;
708  os << "There seems to be something wrong with the quantum numbers of at least one line in your *abs_lines*. " <<
709  "Make sure this is a Zeeman line.\nThe upper quantum numbers are: " <<
710  abs_lines_per_species[II][ii].QuantumNumbers().Upper() <<
711  "\nThe lower quantum numbers are: " <<
712  abs_lines_per_species[II][ii].QuantumNumbers().Lower() <<
713  "\nThe entire line information: " << abs_lines_per_species[II][ii];
714  throw std::runtime_error(os.str());
715  }
716  }
717 
718  if (abs(RS_sum-1.)>margin) //Reasonable confidence?
719  {
720  ostringstream os;
721  os << "The sum of relative strengths is not close to one. This is severly problematic and "
722  "you should look into why this happens.\nIt is currently " << RS_sum << " with DJ: "<<DJ<<", DMain: "<<DMain<<" for line: "<<
723  abs_lines_per_species[II][ii] <<"\n";
724  throw std::runtime_error(os.str());
725  }
726 
727  if(abs_species[II][0].LineMixingType() != SpeciesTag::LINE_MIXING_TYPE_NONE)
728  {
729  for(Index n=0;n<number_M_sm;n++)
730  temp_lut_sm.push_back(line_mixing_data_lut[II][ii]);
731  for(Index n=0;n<number_M_pi;n++)
732  temp_lut_pi.push_back(line_mixing_data_lut[II][ii]);
733  for(Index n=0;n<number_M_sp;n++)
734  temp_lut_sp.push_back(line_mixing_data_lut[II][ii]);
735  }
736 
737  }
738  else
739  {
740  ostringstream os;
741  os << "There are undefined quantum numbers in the line: " << abs_lines_per_species[II][ii] << "\nJ is "<<J<<" and Main is "<<Main<<std::endl;
742  throw std::runtime_error(os.str());
743  }
744  }
745  Tensor3 part_abs_mat(f_grid.nelem(), 4, 4);
746 
747  // Add Pi contribution to final propmat_clearsky
748  xsec_species_line_mixing_wrapper_with_zeeman( part_abs_mat, abs_species, abs_lineshape,
749  temp_abs_lines_pi, isotopologue_ratios,
750  abs_vmrs, abs_p, abs_t, f_grid,
751  theta, eta, 0, II,
752  line_mixing_data,
753  line_mixing_data_lut,
754  temp_lut_pi,
755  verbosity );
756  propmat_clearsky(II, joker, joker, joker) += part_abs_mat;
757 
758  // Add Sigma minus contribution to final propmat_clearsky
759  xsec_species_line_mixing_wrapper_with_zeeman( part_abs_mat, abs_species, abs_lineshape,
760  temp_abs_lines_sm, isotopologue_ratios,
761  abs_vmrs, abs_p, abs_t, f_grid,
762  theta, eta, -1, II,
763  line_mixing_data,
764  line_mixing_data_lut,
765  temp_lut_sm,
766  verbosity );
767  propmat_clearsky(II, joker, joker, joker) += part_abs_mat;
768 
769  // Add Sigma plus contribution to final propmat_clearsky
770  xsec_species_line_mixing_wrapper_with_zeeman( part_abs_mat, abs_species, abs_lineshape,
771  temp_abs_lines_sp, isotopologue_ratios,
772  abs_vmrs, abs_p, abs_t, f_grid,
773  theta, eta, 1, II,
774  line_mixing_data,
775  line_mixing_data_lut,
776  temp_lut_sp,
777  verbosity );
778  propmat_clearsky(II, joker, joker, joker) += part_abs_mat;
779 
780  }
781  }
782  else // if the magnetic field is ignored
783  {
784  for(Index II = 0; II<abs_species.nelem(); II++)
785  {
786  // If the species isn't Zeeman, look at the next species.
787  if(!is_zeeman(abs_species[II])) continue;
788 
789  Tensor3 part_abs_mat(f_grid.nelem(), 4, 4);
790 
791  if(abs_species[II][0].LineMixingType() == SpeciesTag::LINE_MIXING_TYPE_NONE){
792  ArrayOfIndex empty;
793  xsec_species_line_mixing_wrapper_with_zeeman( part_abs_mat, abs_species, abs_lineshape,
794  abs_lines_per_species[II], isotopologue_ratios,
795  abs_vmrs, abs_p, abs_t, f_grid,
796  0,0,1023, II,
797  line_mixing_data,
798  line_mixing_data_lut,
799  empty,
800  verbosity );
801  }
802  else{
803  xsec_species_line_mixing_wrapper_with_zeeman( part_abs_mat, abs_species, abs_lineshape,
804  abs_lines_per_species[II], isotopologue_ratios,
805  abs_vmrs, abs_p, abs_t, f_grid,
806  0,0,1023, II,
807  line_mixing_data,
808  line_mixing_data_lut,
809  line_mixing_data_lut[II],
810  verbosity );
811  }
812  propmat_clearsky(II, joker, joker, joker) += part_abs_mat;
813  }
814  }
815 }
Matrix
The Matrix class.
Definition: matpackI.h:788
frequency_change
Numeric(* frequency_change)(const Rational &, const Rational &, const Rational &, const Numeric &, const Index &, const Index &, const Index &, const Numeric &, const Numeric &)
Definition: m_zeeman.cc:239
QN_J
@ QN_J
Definition: quantum.h:42
MatrixView
The MatrixView class.
Definition: matpackI.h:679
auto_md.h
rational.h
Contains the rational class definition.
absorption.h
Declarations required for the calculation of absorption coefficients.
Tensor3
The Tensor3 class.
Definition: matpackIII.h:348
joker
const Joker joker
SpeciesTag::LINE_MIXING_TYPE_NONE
@ LINE_MIXING_TYPE_NONE
Definition: abs_species_tags.h:130
QuantumNumberRecord::Upper
Rational Upper(Index i) const
Get upper quantum number.
Definition: quantum.h:106
proj
void proj(Vector &c, ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1787
frequency_change_casea
Numeric frequency_change_casea(const Rational &omega, const Rational &m, const Rational &j, const Numeric &Sigma, const Index &DJ, const Index &DM, const Index &Domega, const Numeric &H_mag, const Numeric &GS)
Definition: m_zeeman.cc:297
propmat_clearskyAddZeeman
void propmat_clearskyAddZeeman(Tensor4 &propmat_clearsky, const Vector &f_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfArrayOfLineRecord &abs_lines_per_species, const ArrayOfLineshapeSpec &abs_lineshape, const SpeciesAuxData &isotopologue_ratios, const SpeciesAuxData &isotopologue_quantum, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Vector &rtp_vmr, const Vector &rtp_mag, const Vector &ppath_los, const Index &atmosphere_dim, const ArrayOfArrayOfLineMixingRecord &line_mixing_data, const ArrayOfArrayOfIndex &line_mixing_data_lut, const Index &manual_zeeman_tag, const Numeric &manual_zeeman_magnetic_field_strength, const Numeric &manual_zeeman_theta, const Numeric &manual_zeeman_eta, const Verbosity &verbosity)
WORKSPACE METHOD: propmat_clearskyAddZeeman.
Definition: m_zeeman.cc:365
PI
const Numeric PI
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
QN_Omega
@ QN_Omega
Definition: quantum.h:46
Tensor4
The Tensor4 class.
Definition: matpackIV.h:383
SpeciesAuxData
Auxiliary data for isotopologues.
Definition: absorption.h:440
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:143
dx
#define dx
Definition: continua.cc:21928
cross3
void cross3(VectorView c, const ConstVectorView &a, const ConstVectorView &b)
cross3
Definition: matpackI.cc:1743
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:232
Array
This can be used to make arrays out of anything.
Definition: array.h:107
number_density
Numeric number_density(const Numeric &p, const Numeric &t)
number_density
Definition: physics_funcs.cc:237
is_zeeman
bool is_zeeman(const ArrayOfSpeciesTag &tg)
Is this a Zeeman tag group?
Definition: abs_species_tags.cc:775
LineRecord::QuantumNumbers
const QuantumNumberRecord & QuantumNumbers() const
Quantum numbers.
Definition: linerecord.h:538
checkIsotopologueRatios
void checkIsotopologueRatios(const ArrayOfArrayOfSpeciesTag &abs_species, const SpeciesAuxData &isoratios)
Check that isotopologue ratios for the given species are correctly defined.
Definition: absorption.cc:245
mult
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
CREATE_OUT3
#define CREATE_OUT3
Definition: messages.h:214
SpeciesAuxData::getParam
Numeric getParam(Index species, Index isotopologue, Index col) const
Get single parameter value.
Definition: absorption.h:449
messages.h
Declarations having to do with the four output streams.
QuantumNumberRecord::Lower
Rational Lower(Index i) const
Get lower quantum number.
Definition: quantum.h:103
Rational::toNumeric
Numeric toNumeric() const
Definition: rational.h:52
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
frequency_change_caseb
Numeric frequency_change_caseb(const Rational &n, const Rational &m, const Rational &j, const Numeric &S, const Index &DJ, const Index &DM, const Index &DN, const Numeric &H_mag, const Numeric &GS)
Definition: m_zeeman.cc:261
abs
#define abs(x)
Definition: continua.cc:20458
relative_strength
Numeric relative_strength(const Rational &m, const Rational &j, const Index &dj, const Index &dm)
Definition: m_zeeman.cc:169
LineRecord::Isotopologue
Index Isotopologue() const
The index of the isotopologue species that this line belongs to.
Definition: linerecord.h:307
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
physics_funcs.h
QN_N
@ QN_N
Definition: quantum.h:43
ConstTensor4View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:81
PLANCK_CONST
const Numeric PLANCK_CONST
LANDE_GS
const Numeric LANDE_GS
AbsInputFromRteScalars
void AbsInputFromRteScalars(Vector &abs_p, Vector &abs_t, Matrix &abs_vmrs, const Numeric &rtp_pressure, const Numeric &rtp_temperature, const Vector &rtp_vmr, const Verbosity &)
WORKSPACE METHOD: AbsInputFromRteScalars.
Definition: m_abs.cc:73
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
Verbosity
Definition: messages.h:50
gs_caseb
Numeric gs_caseb(const Numeric &N, const Numeric &J, const Numeric &S, const Numeric &GS)
Definition: m_zeeman.cc:151
ConstTensor4View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:69
matpackIII.h
ConstTensor4View::nbooks
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:63
math_funcs.h
LineRecord
Spectral line catalog data.
Definition: linerecord.h:196
zaaa2cart
void zaaa2cart(Numeric &dx, Numeric &dy, Numeric &dz, const Numeric &za, const Numeric &aa)
zaaa2cart
Definition: ppath.cc:484
LineRecord::setF
void setF(Numeric new_mf)
Set the line center frequency in Hz.
Definition: linerecord.h:342
LineRecord::Species
Index Species() const
The index of the molecular species that this line belongs to.
Definition: linerecord.h:302
abs_species_tags.h
Header file for stuff related to absorption species tags.
xsec_species_line_mixing_wrapper_with_zeeman
void xsec_species_line_mixing_wrapper_with_zeeman(Tensor3View part_abs_mat, const ArrayOfArrayOfSpeciesTag &abs_species, const ArrayOfLineshapeSpec &abs_lineshape, const ArrayOfLineRecord &lr, const SpeciesAuxData &isotopologue_ratios, const Matrix &abs_vmrs, const Vector &abs_p, const Vector &abs_t, const Vector &f_grid, const Numeric &theta, const Numeric &eta, const Index &DM, const Index &this_species, const ArrayOfArrayOfLineMixingRecord &line_mixing_data, const ArrayOfArrayOfIndex &line_mixing_data_lut, const ArrayOfIndex &temp_species_lut, const Verbosity &verbosity)
Definition: m_zeeman.cc:316
ConstTensor3View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:146
xsec_species_line_mixing_wrapper
void xsec_species_line_mixing_wrapper(MatrixView xsec_attenuation, MatrixView xsec_phase, const ArrayOfArrayOfLineMixingRecord &line_mixing_data, const ArrayOfArrayOfIndex &line_mixing_data_lut, ConstVectorView f_grid, ConstVectorView abs_p, ConstVectorView abs_t, ConstMatrixView all_vmrs, const ArrayOfArrayOfSpeciesTag &abs_species, const Index this_species, const ArrayOfLineRecord &abs_lines, const Index ind_ls, const Index ind_lsn, const Numeric cutoff, const SpeciesAuxData &isotopologue_ratios, const Verbosity &verbosity)
This will work as a wrapper for linemixing when abs_species contain relevant data.
Definition: absorption.cc:1530
ConstTensor4View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:75
N
#define N
Definition: rng.cc:195
RAD2DEG
const Numeric RAD2DEG
Rational::isUndefined
bool isUndefined() const
Definition: rational.h:48
ppath.h
Propagation path structure and functions.
mirror_los
void mirror_los(Vector &los_mirrored, ConstVectorView los, const Index &atmosphere_dim)
mirror_los
Definition: rte.cc:2389
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:149
gs_casea
Numeric gs_casea(const Numeric &Omega, const Numeric &J, const Numeric &Sigma, const Numeric &GS)
Definition: m_zeeman.cc:152
rte.h
Declaration of functions in rte.cc.
vector_angle
Numeric vector_angle(ConstVectorView a, ConstVectorView b)
Definition: matpackI.cc:1763
attenuation_matrix
void attenuation_matrix(MatrixView K, const Numeric &theta, const Numeric &eta, const Index &DM)
Definition: m_zeeman.cc:108
DEG2RAD
const Numeric DEG2RAD
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
BOHR_MAGNETON
const Numeric BOHR_MAGNETON
M
#define M
Definition: rng.cc:196
Vector
The Vector class.
Definition: matpackI.h:556
phase_matrix
void phase_matrix(MatrixView K, const Numeric &theta, const Numeric &eta, const Index &DM)
Definition: m_zeeman.cc:55
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:176
LineRecord::setI0
void setI0(Numeric new_mi0)
Set Intensity.
Definition: linerecord.h:365
Rational
Definition: rational.h:35