74 Index& doit_za_grid_size,
78 const Index& N_za_grid,
79 const Index& N_aa_grid,
80 const String& za_grid_opt_file,
89 throw runtime_error(
"N_za_grid must be greater than 15 for accurate results");
90 else if (N_za_grid > 100)
93 out1 <<
"Warning: N_za_grid is very large which means that the \n"
94 <<
"calculation will be very slow.\n";
98 throw runtime_error(
"N_aa_grid must be greater than 5 for accurate results");
99 else if (N_aa_grid > 100)
102 out1 <<
"Warning: N_aa_grid is very large which means that the \n"
103 <<
"calculation will be very slow.\n";
108 nlinspace(scat_aa_grid, 0, 360, N_aa_grid);
112 doit_za_grid_size = N_za_grid;
114 if( za_grid_opt_file ==
"" )
115 nlinspace(scat_za_grid, 0, 180, N_za_grid);
124 Index& doit_conv_flag,
125 Index& doit_iteration_counter,
128 const Tensor6& doit_i_field_old,
131 const Index& max_iterations,
132 const Index& throw_nonconv_error,
139 if( doit_conv_flag != 0 )
140 throw runtime_error(
"Convergence flag is non-zero, which means that this\n"
141 "WSM is not used correctly. *doit_conv_flagAbs* should\n"
142 "be used only in *doit_conv_test_agenda*\n");
149 const Index stokes_dim = doit_i_field.
ncols();
152 if ( epsilon.
nelem() != stokes_dim )
154 "You have to specify limiting values for the "
155 "convergence test for each Stokes component "
156 "separately. That means that *epsilon* must "
157 "have *stokes_dim* elements!"
162 N_p, N_lat, N_lon, N_za, N_aa, stokes_dim))
163 throw runtime_error(
"The fields (Tensor6) *doit_i_field* and \n"
164 "*doit_i_field_old* which are compared in the \n"
165 "convergence test do not have the same size.\n");
170 doit_iteration_counter +=1;
171 out2 <<
" Number of DOIT iteration: " << doit_iteration_counter <<
"\n";
173 if (doit_iteration_counter > max_iterations)
176 out <<
"Method does not converge (number of iterations \n"
177 <<
"is > " << max_iterations <<
"). Either the cloud "
178 <<
"particle number density \n"
179 <<
"is too large or the numerical setup for the DOIT \n"
180 <<
"calculation is not correct. In case of limb \n"
181 <<
"simulations please make sure that you use an \n"
182 <<
"optimized zenith angle grid. \n"
183 <<
"*doit_i_field* might be wrong.\n";
184 if( throw_nonconv_error != 0)
191 out1 <<
"Warning in DOIT calculation (output set to NaN):\n"
198 out1 <<
"Warning in DOIT calculation (output equals current status):\n"
205 for (
Index p_index = 0; p_index < N_p; p_index++)
207 for (
Index lat_index = 0; lat_index < N_lat; lat_index++)
209 for (
Index lon_index = 0; lon_index <N_lon; lon_index++)
211 for (
Index scat_za_index = 0; scat_za_index < N_za;
214 for (
Index scat_aa_index = 0; scat_aa_index < N_aa;
217 for (
Index stokes_index = 0; stokes_index <
218 stokes_dim; stokes_index ++)
221 (doit_i_field(p_index, lat_index, lon_index,
222 scat_za_index, scat_aa_index,
224 doit_i_field_old(p_index, lat_index,
225 lon_index, scat_za_index,
233 if(
abs(diff) > epsilon[stokes_index])
235 out1 <<
"difference: " << diff <<
"\n";
255 Index& doit_conv_flag,
256 Index& doit_iteration_counter,
259 const Tensor6& doit_i_field_old,
261 const Index& f_index,
264 const Index& max_iterations,
265 const Index& throw_nonconv_error,
273 if( doit_conv_flag != 0 )
274 throw runtime_error(
"Convergence flag is non-zero, which means that this \n"
275 "WSM is not used correctly. *doit_conv_flagAbs* should\n"
276 "be used only in *doit_conv_test_agenda*\n");
283 const Index stokes_dim = doit_i_field.
ncols();
286 if ( epsilon.
nelem() != stokes_dim )
288 "You have to specify limiting values for the "
289 "convergence test for each Stokes component "
290 "separately. That means that *epsilon* must "
291 "have *stokes_dim* elements!"
296 N_p, N_lat, N_lon, N_za, N_aa, stokes_dim))
297 throw runtime_error(
"The fields (Tensor6) *doit_i_field* and \n"
298 "*doit_i_field_old* which are compared in the \n"
299 "convergence test do not have the same size.\n");
303 if( f_grid.
nelem() == 0 )
304 throw runtime_error(
"The frequency grid is empty." );
308 if ( f_index >= f_grid.
nelem() )
309 throw runtime_error(
"*f_index* is greater than number of elements in the\n"
310 "frequency grid.\n");
314 doit_iteration_counter +=1;
315 out2 <<
" Number of DOIT iteration: " << doit_iteration_counter <<
"\n";
317 if (doit_iteration_counter > max_iterations)
320 out <<
"At frequency " << f_grid[f_index] <<
" GHz \n"
321 <<
"method does not converge (number of iterations \n"
322 <<
"is > " << max_iterations <<
"). Either the cloud particle"
323 <<
" number density \n"
324 <<
"is too large or the numerical setup for the DOIT \n"
325 <<
"calculation is not correct. In case of limb \n"
326 <<
"simulations please make sure that you use an \n"
327 <<
"optimized zenith angle grid. \n"
328 <<
"*doit_i_field* might be wrong.\n";
329 if( throw_nonconv_error != 0)
336 out1 <<
"Warning in DOIT calculation (output set to NaN):\n"
343 out1 <<
"Warning in DOIT calculation (output equals current status):\n"
350 for (
Index p_index = 0; p_index < N_p; p_index++)
352 for (
Index lat_index = 0; lat_index < N_lat; lat_index++)
354 for (
Index lon_index = 0; lon_index <N_lon; lon_index++)
356 for (
Index scat_za_index = 0; scat_za_index < N_za;
359 for (
Index scat_aa_index = 0; scat_aa_index < N_aa;
362 for (
Index stokes_index = 0; stokes_index <
363 stokes_dim; stokes_index ++)
366 doit_i_field(p_index, lat_index, lon_index,
367 scat_za_index, scat_aa_index,
369 doit_i_field_old(p_index, lat_index,
370 lon_index, scat_za_index,
378 if(
abs(diff_bt) > epsilon[stokes_index])
380 out1 <<
"BT difference: " << diff_bt
381 <<
" in stokes dim " << stokes_index <<
"\n";
399 Index& doit_conv_flag,
400 Index& doit_iteration_counter,
403 const Tensor6& doit_i_field_old,
405 const Index& f_index,
408 const Index& max_iterations,
409 const Index& throw_nonconv_error,
417 if( doit_conv_flag != 0 )
418 throw runtime_error(
"Convergence flag is non-zero, which means that this \n"
419 "WSM is not used correctly. *doit_conv_flagAbs* should\n"
420 "be used only in *doit_conv_test_agenda*\n");
427 const Index stokes_dim = doit_i_field.
ncols();
430 if ( epsilon.
nelem() != stokes_dim )
432 "You have to specify limiting values for the "
433 "convergence test for each Stokes component "
434 "separately. That means that *epsilon* must "
435 "have *stokes_dim* elements!"
440 N_p, N_lat, N_lon, N_za, N_aa, stokes_dim))
441 throw runtime_error(
"The fields (Tensor6) *doit_i_field* and \n"
442 "*doit_i_field_old* which are compared in the \n"
443 "convergence test do not have the same size.\n");
447 if( f_grid.
nelem() == 0 )
448 throw runtime_error(
"The frequency grid is empty." );
452 if ( f_index >= f_grid.
nelem() )
453 throw runtime_error(
"*f_index* is greater than number of elements in the\n"
454 "frequency grid.\n");
459 doit_iteration_counter +=1;
460 out2 <<
" Number of DOIT iteration: " << doit_iteration_counter <<
"\n";
462 if (doit_iteration_counter > max_iterations)
465 out <<
"Method does not converge (number of iterations \n"
466 <<
"is > " << max_iterations <<
"). Either the cloud"
467 <<
" particle number density \n"
468 <<
"is too large or the numerical setup for the DOIT \n"
469 <<
"calculation is not correct. In case of limb \n"
470 <<
"simulations please make sure that you use an \n"
471 <<
"optimized zenith angle grid. \n";
472 if( throw_nonconv_error != 0)
479 out1 <<
"Warning in DOIT calculation (output set to NaN):\n"
486 out1 <<
"Warning in DOIT calculation (output equals current status):\n"
499 for (
Index p_index = 0; p_index < N_p; p_index++)
501 for (
Index lat_index = 0; lat_index < N_lat; lat_index++)
503 for (
Index lon_index = 0; lon_index <N_lon; lon_index++)
505 for (
Index scat_za_index = 0; scat_za_index < N_za;
508 for (
Index scat_aa_index = 0; scat_aa_index < N_aa;
513 doit_i_field(p_index, lat_index,
515 scat_za_index, scat_aa_index, i) -
516 doit_i_field_old(p_index, lat_index,
517 lon_index, scat_za_index,
526 lqs[i] = sqrt(lqs[i]);
527 lqs[i] /= (
Numeric)(N_p*N_lat*N_lon*N_za*N_aa);
532 if (lqs[i] >= epsilon[i] )
536 out1 <<
"lqs [I]: " << lqs[0] <<
"\n";
546 const Agenda& doit_scat_field_agenda,
547 const Agenda& doit_rte_agenda,
548 const Agenda& doit_conv_test_agenda,
554 chk_not_empty(
"doit_scat_field_agenda", doit_scat_field_agenda);
556 chk_not_empty(
"doit_conv_test_agenda", doit_conv_test_agenda);
563 Tensor6 doit_i_field_old_local;
564 Index doit_conv_flag_local;
565 Index doit_iteration_counter_local;
572 doit_i_field.
nrows(), doit_i_field.
ncols(), 0.);
574 doit_conv_flag_local = 0;
575 doit_iteration_counter_local = 0;
577 while(doit_conv_flag_local == 0) {
580 doit_i_field_old_local = doit_i_field;
585 out2 <<
" Execute doit_scat_field_agenda. \n";
588 doit_scat_field_agenda);
591 out2 <<
" Execute doit_rte_agenda. \n";
597 doit_iteration_counter_local,
599 doit_i_field_old_local,
600 doit_conv_test_agenda);
612 const Tensor6& doit_scat_field,
615 const Agenda& propmat_clearsky_agenda,
618 const Agenda& spt_calc_agenda,
619 const Vector& scat_za_grid,
622 const Agenda& opt_prop_part_agenda,
624 const Agenda& ppath_step_agenda,
625 const Numeric& ppath_lraytrace,
628 const Vector& refellipsoid,
632 const Index& f_index,
633 const Agenda& surface_rtprop_agenda,
634 const Index& doit_za_interp,
641 out2 <<
" doit_i_fieldUpdate1D: Radiative transfer calculation in cloudbox\n";
642 out2 <<
" ------------------------------------------------------------- \n";
648 chk_not_empty(
"opt_prop_part_agenda", opt_prop_part_agenda);
651 if (cloudbox_limits.
nelem() != 2)
653 "The cloudbox dimension is not 1D! \n"
654 "Do you really want to do a 1D calculation? \n"
655 "If not, use *doit_i_fieldUpdateSeq3D*.\n"
661 if (scat_za_grid[0] != 0. || scat_za_grid[N_scat_za-1] != 180.)
662 throw runtime_error(
"The range of *scat_za_grid* must [0 180].");
664 if( p_grid.
nelem() < 2 )
665 throw runtime_error(
"The length of *p_grid* must be >= 2." );
674 if( f_grid.
nelem() == 0 )
675 throw runtime_error(
"The frequency grid is empty." );
679 if ( f_index >= f_grid.
nelem() )
680 throw runtime_error(
"*f_index* is greater than number of elements in the\n"
681 "frequency grid.\n");
683 if( !(doit_za_interp == 0 || doit_za_interp == 1 ) )
684 throw runtime_error(
"Interpolation method is not defined. Use \n"
685 "*doit_za_interpSet*.\n");
687 const Index stokes_dim = doit_scat_field.
ncols();
688 assert(stokes_dim > 0 || stokes_dim < 5);
693 (cloudbox_limits[1] - cloudbox_limits[0]) + 1, 1, 1,
694 N_scat_za, 1, stokes_dim));
696 assert(
is_size( doit_scat_field,
697 (cloudbox_limits[1] - cloudbox_limits[0]) + 1, 1, 1,
698 N_scat_za, 1, stokes_dim));
708 out3 <<
"Calculate single particle properties \n";
718 Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1,
719 stokes_dim, stokes_dim, 0.);
720 Tensor4 abs_vec_field(cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1,
723 Tensor6 doit_i_field_old(doit_i_field);
726 Index scat_aa_index_local = 0;
729 for(
Index scat_za_index_local = 0; scat_za_index_local < N_scat_za;
730 scat_za_index_local ++)
739 opt_prop_part_agenda, scat_za_index_local,
741 cloudbox_limits, t_field, pnd_field, verbosity);
747 for(
Index p_index = cloudbox_limits[0]; p_index
748 <= cloudbox_limits[1]; p_index ++)
750 if ( (p_index!=0) || (scat_za_grid[scat_za_index_local] <= 90.))
753 p_index, scat_za_index_local,
755 cloudbox_limits, doit_i_field_old,
757 propmat_clearsky_agenda, vmr_field,
758 ppath_step_agenda, ppath_lraytrace,
759 p_grid, z_field, refellipsoid,
760 t_field, f_grid, f_index, ext_mat_field,
762 surface_rtprop_agenda, doit_za_interp,
779 const Agenda& propmat_clearsky_agenda,
782 const Agenda& spt_calc_agenda,
783 const Vector& scat_za_grid,
784 const Vector& scat_aa_grid,
787 const Agenda& opt_prop_part_agenda,
789 const Agenda& ppath_step_agenda,
790 const Numeric& ppath_lraytrace,
793 const Vector& refellipsoid,
797 const Index& f_index,
798 const Agenda& surface_rtprop_agenda,
799 const Index& doit_za_interp,
800 const Index& normalize,
801 const Numeric& norm_error_threshold,
802 const Index& norm_debug,
808 out2<<
" doit_i_fieldUpdateSeq1D: Radiative transfer calculation in cloudbox\n";
809 out2 <<
" ------------------------------------------------------------- \n";
814 chk_not_empty(
"propmat_clearsky_agenda", propmat_clearsky_agenda);
816 chk_not_empty(
"opt_prop_part_agenda", opt_prop_part_agenda);
819 if (cloudbox_limits.
nelem() != 2)
821 "The cloudbox dimension is not 1D! \n"
822 "Do you really want to do a 1D calculation? \n"
823 "For 3D use *doit_i_fieldUpdateSeq3D*.\n"
829 if (scat_za_grid[0] != 0. || scat_za_grid[N_scat_za-1] != 180.)
830 throw runtime_error(
"The range of *scat_za_grid* must [0 180].");
832 if( p_grid.
nelem() < 2 )
833 throw runtime_error(
"The length of *p_grid* must be >= 2." );
842 if( f_grid.
nelem() == 0 )
843 throw runtime_error(
"The frequency grid is empty." );
847 if ( f_index >= f_grid.
nelem() )
848 throw runtime_error(
"*f_index* is greater than number of elements in the\n"
849 "frequency grid.\n");
851 if( !(doit_za_interp == 0 || doit_za_interp == 1 ) )
852 throw runtime_error(
"Interpolation method is not defined. Use \n"
853 "*doit_za_interpSet*.\n");
855 const Index stokes_dim = doit_scat_field.
ncols();
856 assert(stokes_dim > 0 || stokes_dim < 5);
861 (cloudbox_limits[1] - cloudbox_limits[0]) + 1, 1, 1,
862 N_scat_za, 1, stokes_dim));
864 assert(
is_size( doit_scat_field,
865 (cloudbox_limits[1] - cloudbox_limits[0]) + 1, 1, 1,
866 N_scat_za, 1, stokes_dim));
876 out3 <<
"Calculate single particle properties \n";
886 Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1,
887 stokes_dim, stokes_dim, 0.);
888 Tensor4 abs_vec_field(cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1,
894 Numeric theta_lim = 180. - asin((refellipsoid[0]+
895 z_field(cloudbox_limits[0],0,0))/
897 z_field(cloudbox_limits[1],0,0)))*
RAD2DEG;
909 Index scat_aa_index_local = 0;
920 scat_za_grid, scat_aa_grid,
922 opt_prop_part_agenda,
924 norm_error_threshold,
930 for(
Index scat_za_index_local = 0; scat_za_index_local < N_scat_za;
931 scat_za_index_local ++)
939 spt_calc_agenda, opt_prop_part_agenda,
940 scat_za_index_local, scat_aa_index_local,
941 cloudbox_limits, t_field, pnd_field, verbosity);
950 if ( scat_za_grid[scat_za_index_local] <= 90.)
957 for(
Index p_index = cloudbox_limits[1]-1; p_index
958 >= cloudbox_limits[0]; p_index --)
961 p_index, scat_za_index_local, scat_za_grid,
962 cloudbox_limits, doit_scat_field,
963 propmat_clearsky_agenda, vmr_field,
964 ppath_step_agenda, ppath_lraytrace,
965 p_grid, z_field, refellipsoid,
966 t_field, f_grid, f_index,
967 ext_mat_field, abs_vec_field,
968 surface_rtprop_agenda, doit_za_interp,
972 else if ( scat_za_grid[scat_za_index_local] >= theta_lim)
977 for(
Index p_index = cloudbox_limits[0]+1; p_index
978 <= cloudbox_limits[1]; p_index ++)
981 p_index, scat_za_index_local, scat_za_grid,
982 cloudbox_limits, doit_scat_field,
983 propmat_clearsky_agenda, vmr_field,
984 ppath_step_agenda, ppath_lraytrace,
985 p_grid, z_field, refellipsoid,
986 t_field, f_grid, f_index,
987 ext_mat_field, abs_vec_field,
988 surface_rtprop_agenda, doit_za_interp,
1002 bool conv_flag =
false;
1004 while (!conv_flag && limb_it < 10)
1007 doit_i_field_limb = doit_i_field(
joker, 0, 0, scat_za_index_local, 0,
joker);
1008 for(
Index p_index = cloudbox_limits[0];
1009 p_index <= cloudbox_limits[1]; p_index ++)
1018 p_index, scat_za_index_local,
1020 cloudbox_limits, doit_scat_field,
1021 propmat_clearsky_agenda, vmr_field,
1022 ppath_step_agenda, ppath_lraytrace,
1023 p_grid, z_field, refellipsoid,
1024 t_field, f_grid, f_index,
1025 ext_mat_field, abs_vec_field,
1026 surface_rtprop_agenda, doit_za_interp,
1033 for (
Index p_index = 0;
1034 conv_flag && p_index < doit_i_field.
nvitrines(); p_index++)
1036 for (
Index stokes_index = 0;
1037 conv_flag && stokes_index < stokes_dim; stokes_index ++)
1040 doit_i_field(p_index, 0, 0, scat_za_index_local, 0, stokes_index)
1041 - doit_i_field_limb(p_index, stokes_index);
1047 if (
abs(diff_bt) > epsilon[stokes_index])
1049 out2 <<
"Limb BT difference: " << diff_bt
1050 <<
" in stokes dim " << stokes_index <<
"\n";
1056 out2 <<
"Limb iterations: " << limb_it <<
"\n";
1068 const Tensor6& doit_scat_field,
1071 const Agenda& propmat_clearsky_agenda,
1074 const Agenda& spt_calc_agenda,
1075 const Vector& scat_za_grid,
1076 const Vector& scat_aa_grid,
1079 const Agenda& opt_prop_part_agenda,
1081 const Agenda& ppath_step_agenda,
1082 const Numeric& ppath_lraytrace,
1087 const Vector& refellipsoid,
1091 const Index& f_index,
1092 const Index& doit_za_interp,
1098 out2<<
" doit_i_fieldUpdateSeq3D: Radiative transfer calculatiuon in cloudbox.\n";
1099 out2 <<
" ------------------------------------------------------------- \n";
1104 chk_not_empty(
"propmat_clearsky_agenda",propmat_clearsky_agenda);
1106 chk_not_empty(
"opt_prop_part_agenda", opt_prop_part_agenda);
1109 if (cloudbox_limits.
nelem() != 6)
1110 throw runtime_error(
1111 "The cloudbox dimension is not 3D! \n"
1112 "Do you really want to do a 3D calculation? \n"
1113 "For 1D use *doit_i_fieldUpdateSeq1D*.\n"
1117 const Index N_scat_za = scat_za_grid.
nelem();
1119 if (scat_za_grid[0] != 0. || scat_za_grid[N_scat_za-1] != 180.)
1120 throw runtime_error(
"The range of *scat_za_grid* must [0 180].");
1123 const Index N_scat_aa = scat_aa_grid.
nelem();
1125 if (scat_aa_grid[0] != 0. || scat_aa_grid[N_scat_aa-1] != 360.)
1126 throw runtime_error(
"The range of *scat_za_grid* must [0 360].");
1139 if( f_grid.
nelem() == 0 )
1140 throw runtime_error(
"The frequency grid is empty." );
1144 if ( f_index >= f_grid.
nelem() )
1145 throw runtime_error(
"*f_index* is greater than number of elements in the\n"
1146 "frequency grid.\n");
1148 if( !(doit_za_interp == 0 || doit_za_interp == 1 ) )
1149 throw runtime_error(
"Interpolation method is not defined. Use \n"
1150 "*doit_za_interpSet*.\n");
1152 const Index stokes_dim = doit_scat_field.
ncols();
1153 assert(stokes_dim > 0 || stokes_dim < 5);
1156 assert(
is_size( doit_i_field,
1157 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1158 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
1159 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
1164 assert(
is_size( doit_scat_field,
1165 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1166 (cloudbox_limits[3] - cloudbox_limits[2]) + 1,
1167 (cloudbox_limits[5] - cloudbox_limits[4]) + 1,
1180 out3 <<
"Calculate single particle properties \n";
1190 const Index p_low = cloudbox_limits[0];
1191 const Index p_up = cloudbox_limits[1];
1192 const Index lat_low = cloudbox_limits[2];
1193 const Index lat_up = cloudbox_limits[3];
1194 const Index lon_low = cloudbox_limits[4];
1195 const Index lon_up = cloudbox_limits[5];
1199 Tensor5 ext_mat_field(p_up-p_low+1, lat_up-lat_low+1, lon_up-lon_low+1,
1200 stokes_dim, stokes_dim, 0.);
1201 Tensor4 abs_vec_field(p_up-p_low+1, lat_up-lat_low+1, lon_up-lon_low+1,
1206 for(
Index scat_za_index = 0; scat_za_index < N_scat_za; scat_za_index ++)
1210 for(
Index scat_aa_index = 1; scat_aa_index < N_scat_aa; scat_aa_index ++)
1221 opt_prop_part_agenda, scat_za_index,
1222 scat_aa_index, cloudbox_limits, t_field,
1223 pnd_field, verbosity);
1226 Vector stokes_vec(stokes_dim,0.);
1228 Numeric theta_lim = 180. - asin((refellipsoid[0]+z_field(p_low,0,0))
1229 /(refellipsoid[0]+z_field(p_up,0,0)))
1233 if ( scat_za_grid[scat_za_index] <= 90.)
1240 for(
Index p_index = p_up-1; p_index >= p_low; p_index --)
1242 for(
Index lat_index = lat_low; lat_index <= lat_up;
1245 for(
Index lon_index = lon_low; lon_index <= lon_up;
1250 lon_index, scat_za_index,
1251 scat_aa_index, scat_za_grid,
1252 scat_aa_grid, cloudbox_limits,
1254 propmat_clearsky_agenda,
1255 vmr_field, ppath_step_agenda,
1256 ppath_lraytrace, p_grid,
1257 lat_grid, lon_grid, z_field,
1258 refellipsoid, t_field,
1260 ext_mat_field, abs_vec_field,
1267 else if ( scat_za_grid[scat_za_index] > theta_lim)
1272 for(
Index p_index = p_low+1; p_index <= p_up; p_index ++)
1274 for(
Index lat_index = lat_low; lat_index <= lat_up;
1277 for(
Index lon_index = lon_low; lon_index <= lon_up;
1282 lon_index, scat_za_index,
1283 scat_aa_index, scat_za_grid,
1284 scat_aa_grid, cloudbox_limits,
1286 propmat_clearsky_agenda,
1287 vmr_field, ppath_step_agenda,
1288 ppath_lraytrace, p_grid,
1289 lat_grid, lon_grid, z_field,
1290 refellipsoid, t_field,
1292 ext_mat_field, abs_vec_field,
1307 else if ( scat_za_grid[scat_za_index] > 90. &&
1308 scat_za_grid[scat_za_index] < theta_lim )
1310 for(
Index p_index = p_low; p_index <= p_up; p_index ++)
1316 if (!(p_index == 0 && scat_za_grid[scat_za_index] > 90.))
1318 for(
Index lat_index = lat_low; lat_index <= lat_up;
1321 for(
Index lon_index = lon_low; lon_index <= lon_up;
1327 lon_index, scat_za_index,
1333 propmat_clearsky_agenda,
1334 vmr_field, ppath_step_agenda,
1335 ppath_lraytrace, p_grid,
1339 t_field, f_grid, f_index,
1364 Index& scat_za_index ,
1366 const Tensor6& doit_scat_field,
1369 const Agenda& propmat_clearsky_agenda,
1372 const Agenda& spt_calc_agenda,
1373 const Vector& scat_za_grid,
1376 const Agenda& opt_prop_part_agenda,
1383 const Index& f_index,
1389 out2 <<
" doit_i_fieldUpdateSeq1DPP: Radiative transfer calculation in cloudbox.\n";
1390 out2 <<
" --------------------------------------------------------------------- \n";
1392 const Index stokes_dim = doit_scat_field.
ncols();
1397 if (stokes_dim < 0 || stokes_dim > 4)
1398 throw runtime_error(
1399 "The dimension of stokes vector must be"
1402 assert(
is_size( doit_i_field,
1403 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1406 scat_za_grid.
nelem(),
1410 assert(
is_size( doit_scat_field,
1411 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
1414 scat_za_grid.
nelem(),
1419 assert( f_index <= f_grid.
nelem() );
1426 const Index N_scat_za = scat_za_grid.
nelem();
1433 out3 <<
"Calculate single particle properties \n";
1444 Tensor5 ext_mat_field(cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1,
1445 stokes_dim, stokes_dim, 0.);
1446 Tensor4 abs_vec_field(cloudbox_limits[1] - cloudbox_limits[0] + 1, 1, 1,
1450 for(scat_za_index = 0; scat_za_index < N_scat_za; scat_za_index ++)
1454 Index scat_aa_index = 0;
1458 opt_prop_part_agenda, scat_za_index, scat_aa_index,
1459 cloudbox_limits, t_field,
1460 pnd_field, verbosity);
1466 Vector stokes_vec(stokes_dim,0.);
1468 if ( scat_za_grid[scat_za_index] <= 90)
1478 for(
Index p_index = cloudbox_limits[1] -1; p_index
1479 >= cloudbox_limits[0]; p_index --)
1482 p_index, scat_za_index,
1486 propmat_clearsky_agenda,
1496 else if ( scat_za_grid[scat_za_index] > 90)
1501 for(
Index p_index = cloudbox_limits[0]+1; p_index
1502 <= cloudbox_limits[1]; p_index ++)
1505 p_index, scat_za_index,
1509 propmat_clearsky_agenda,
1526 Index& scat_p_index,
1527 Index& scat_lat_index,
1528 Index& scat_lon_index,
1529 Index& scat_za_index,
1530 Index& scat_aa_index,
1533 Tensor4& doit_i_field1D_spectrum,
1537 Index& doit_is_initialized,
1539 const Index& stokes_dim,
1540 const Index& atmosphere_dim,
1542 const Vector& scat_za_grid,
1543 const Vector& scat_aa_grid,
1544 const Index& doit_za_grid_size,
1545 const Index& cloudbox_on,
1553 doit_is_initialized = 0;
1554 out0 <<
" Cloudbox is off, DOIT calculation will be skipped.\n";
1560 if (stokes_dim < 0 || stokes_dim > 4)
1561 throw runtime_error(
1562 "The dimension of stokes vector must be"
1568 const Index N_scat_za = scat_za_grid.
nelem();
1570 if (scat_za_grid[0] != 0. || scat_za_grid[N_scat_za-1] != 180.)
1571 throw runtime_error(
"The range of *scat_za_grid* must [0 180].");
1574 throw runtime_error(
"*scat_za_grid* must be increasing.");
1577 const Index N_scat_aa = scat_aa_grid.
nelem();
1579 if (scat_aa_grid[0] != 0. || scat_aa_grid[N_scat_aa-1] != 360.)
1580 throw runtime_error(
"The range of *scat_aa_grid* must [0 360].");
1582 if (doit_za_grid_size < 16)
1583 throw runtime_error(
1584 "*doit_za_grid_size* must be greater than 15 for accurate results");
1585 else if (doit_za_grid_size > 100)
1588 out1 <<
"Warning: doit_za_grid_size is very large which means that the \n"
1589 <<
"calculation will be very slow. The recommended value is 19.\n";
1592 if ( cloudbox_limits.
nelem()!= 2*atmosphere_dim)
1593 throw runtime_error(
1594 "*cloudbox_limits* is a vector which contains the"
1595 "upper and lower limit of the cloud for all "
1596 "atmospheric dimensions. So its dimension must"
1597 "be 2 x *atmosphere_dim*");
1599 if (scat_data_array.
nelem() == 0)
1600 throw runtime_error(
1601 "No scattering data files have been added.\n"
1602 "Please use the WSM *ParticleTypeAdd* or \n"
1603 "*ParticleTypeAddAll* to define the cloud \n"
1604 "properties for the scattering calculation.\n"
1620 if (atmosphere_dim == 1)
1622 doit_i_field.
resize((cloudbox_limits[1] - cloudbox_limits[0]) +1,
1625 scat_za_grid.
nelem(),
1629 doit_scat_field.
resize((cloudbox_limits[1] - cloudbox_limits[0]) +1,
1632 scat_za_grid.
nelem(),
1636 (cloudbox_limits[1] - cloudbox_limits[0]) +1,
1637 scat_za_grid.
nelem(),
1640 scat_za_grid.
nelem(), 1, stokes_dim);
1642 (cloudbox_limits[1]-cloudbox_limits[0])+1, 2, 0,
1643 scat_za_grid.
nelem(), 1, stokes_dim);
1645 (cloudbox_limits[1]-cloudbox_limits[0])+1, 0, 2,
1646 scat_za_grid.
nelem(), 1, stokes_dim);
1648 else if (atmosphere_dim == 3)
1650 doit_i_field.
resize((cloudbox_limits[1] - cloudbox_limits[0]) +1,
1651 (cloudbox_limits[3] - cloudbox_limits[2]) +1,
1652 (cloudbox_limits[5] - cloudbox_limits[4]) +1,
1653 scat_za_grid.
nelem(),
1654 scat_aa_grid.
nelem(),
1657 doit_scat_field.
resize((cloudbox_limits[1] - cloudbox_limits[0]) +1,
1658 (cloudbox_limits[3] - cloudbox_limits[2]) +1,
1659 (cloudbox_limits[5] - cloudbox_limits[4]) +1,
1660 scat_za_grid.
nelem(),
1661 scat_aa_grid.
nelem(),
1666 throw runtime_error(
1667 "Scattering calculations are not possible for a 2D"
1668 "atmosphere. If you want to do scattering calculations"
1669 "*atmosphere_dim* has to be either 1 or 3"
1674 doit_scat_field = 0.;
1675 doit_i_field1D_spectrum = 0.;
1679 doit_is_initialized = 1;
1685 const Index& doit_iteration_counter,
1698 os << doit_iteration_counter;
1701 if( iterations[0] == 0 )
1712 if (doit_iteration_counter == iterations[i])
1726 const Agenda& pha_mat_spt_agenda,
1730 const Index& atmosphere_dim,
1732 const Vector& scat_za_grid,
1733 const Vector& scat_aa_grid,
1734 const Index& doit_za_grid_size,
1749 if (scat_za_grid[0] != 0. || scat_za_grid[Nza-1] != 180.)
1750 throw runtime_error(
"The range of *scat_za_grid* must [0 180].");
1755 if (scat_aa_grid[0] != 0. || scat_aa_grid[Naa-1] != 360.)
1756 throw runtime_error(
"The range of *scat_za_grid* must [0 360].");
1759 const Index stokes_dim = doit_scat_field.
ncols();
1760 assert(stokes_dim > 0 || stokes_dim < 5);
1770 if (atmosphere_dim == 1)
1773 (cloudbox_limits[1] - cloudbox_limits[0]) +1,
1774 1, 1, Nza, 1, stokes_dim));
1775 assert(
is_size(doit_scat_field,
1776 (cloudbox_limits[1] - cloudbox_limits[0]) +1,
1777 1, 1, scat_za_grid.
nelem(), 1, stokes_dim));
1779 else if (atmosphere_dim == 3)
1781 assert (
is_size( doit_i_field,
1782 (cloudbox_limits[1] - cloudbox_limits[0]) +1,
1783 (cloudbox_limits[3] - cloudbox_limits[2]) +1,
1784 (cloudbox_limits[5] - cloudbox_limits[4]) +1,
1785 Nza, Naa, stokes_dim));
1786 assert (
is_size( doit_scat_field,
1787 (cloudbox_limits[1] - cloudbox_limits[0]) +1,
1788 (cloudbox_limits[3] - cloudbox_limits[2]) +1,
1789 (cloudbox_limits[5] - cloudbox_limits[4]) +1,
1790 Nza, Naa, stokes_dim));
1795 os <<
"The atmospheric dimension must be 1D or 3D \n"
1796 <<
"for scattering calculations using the DOIT\n"
1797 <<
"module, but it is not. The value of *atmosphere_dim*\n"
1798 <<
"is " << atmosphere_dim <<
".";
1799 throw runtime_error( os.str() );
1802 if ( cloudbox_limits.
nelem()!= 2*atmosphere_dim)
1803 throw runtime_error(
1804 "*cloudbox_limits* is a vector which contains the"
1805 "upper and lower limit of the cloud for all "
1806 "atmospheric dimensions. So its dimension must"
1807 "be 2 x *atmosphere_dim*");
1811 if (doit_za_grid_size != Nza)
1812 throw runtime_error(
1813 "The zenith angle grids for the computation of\n"
1814 "the scattering integral and the RT part must \n"
1815 "be equal. Check definitions in \n"
1816 "*DoitAngularGridsSet*. The keyword \n"
1817 "'za_grid_opt_file' should be empty. \n"
1823 Tensor4 pha_mat_local(doit_za_grid_size, scat_aa_grid.
nelem(),
1824 stokes_dim, stokes_dim, 0.);
1826 Tensor5 pha_mat_spt_local(pnd_field.
nbooks(), doit_za_grid_size,
1827 scat_aa_grid.
nelem(), stokes_dim, stokes_dim, 0.);
1831 grid_stepsize[0] = 180./(
Numeric)(doit_za_grid_size - 1);
1832 grid_stepsize[1] = 360./(
Numeric)(Naa - 1);
1834 Tensor3 product_field(Nza, Naa, stokes_dim, 0);
1836 out2 <<
" Calculate the scattered field\n";
1838 if ( atmosphere_dim == 1 )
1840 Index scat_aa_index_local = 0;
1844 for (
Index p_index = 0; p_index<=cloudbox_limits[1]-cloudbox_limits[0] ;
1847 Numeric rtp_temperature_local =
1848 t_field(p_index + cloudbox_limits[0], 0, 0);
1850 for (
Index scat_za_index_local = 0;
1851 scat_za_index_local < Nza; scat_za_index_local ++)
1854 Index index_zero = 0;
1857 out3 <<
"Calculate the phase matrix \n";
1859 scat_za_index_local,
1863 scat_aa_index_local,
1864 rtp_temperature_local,
1865 pha_mat_spt_agenda);
1868 pha_matCalc(pha_mat_local, pha_mat_spt_local, pnd_field,
1869 atmosphere_dim, p_index, 0,
1872 out3 <<
"Multiplication of phase matrix with incoming" <<
1879 for (
Index za_in = 0; za_in < Nza; ++ za_in)
1881 for (
Index aa_in = 0; aa_in < Naa; ++ aa_in)
1886 for (
Index i = 0; i < stokes_dim; i++)
1888 for (
Index j = 0; j< stokes_dim; j++)
1890 product_field(za_in, aa_in, i) +=
1891 pha_mat_local(za_in, aa_in, i, j) *
1892 doit_i_field(p_index, 0, 0, za_in, 0, j);
1900 for (
Index i = 0; i < stokes_dim; i++)
1902 doit_scat_field( p_index, 0, 0, scat_za_index_local, 0, i)
1916 else if( atmosphere_dim == 3 )
1922 for (
Index p_index = 0; p_index <=
1923 cloudbox_limits[1] - cloudbox_limits[0];
1926 for (
Index lat_index = 0; lat_index <=
1927 cloudbox_limits[3] - cloudbox_limits[2]; lat_index++)
1929 for (
Index lon_index = 0; lon_index <=
1930 cloudbox_limits[5]-cloudbox_limits[4]; lon_index++)
1932 Numeric rtp_temperature_local =
1933 t_field(p_index + cloudbox_limits[0],
1934 lat_index + cloudbox_limits[2],
1935 lon_index + cloudbox_limits[4]);
1937 for (
Index scat_aa_index_local = 1;
1938 scat_aa_index_local < Naa;
1939 scat_aa_index_local++)
1941 for (
Index scat_za_index_local = 0;
1942 scat_za_index_local < Nza;
1943 scat_za_index_local ++)
1945 out3 <<
"Calculate phase matrix \n";
1947 scat_za_index_local,
1951 scat_aa_index_local,
1952 rtp_temperature_local,
1953 pha_mat_spt_agenda);
1967 for (
Index za_in = 0; za_in < Nza; ++ za_in)
1969 for (
Index aa_in = 0; aa_in < Naa; ++ aa_in)
1973 for (
Index i = 0; i < stokes_dim; i++)
1975 for (
Index j = 0; j< stokes_dim; j++)
1977 product_field(za_in, aa_in, i) +=
1979 (za_in, aa_in, i, j) *
1980 doit_i_field(p_index, lat_index,
1982 scat_za_index_local,
1983 scat_aa_index_local,
1993 for (
Index i = 0; i < stokes_dim; i++)
1995 doit_scat_field( p_index,
1998 scat_za_index_local,
1999 scat_aa_index_local,
2026 const Agenda& pha_mat_spt_agenda,
2030 const Index& atmosphere_dim,
2032 const Vector& scat_za_grid,
2033 const Vector& scat_aa_grid,
2034 const Index& doit_za_grid_size,
2035 const Index& doit_za_interp,
2049 if (scat_za_grid[0] != 0. || scat_za_grid[Nza-1] != 180.)
2050 throw runtime_error(
"The range of *scat_za_grid* must [0 180].");
2055 if (scat_aa_grid[0] != 0. || scat_aa_grid[Naa-1] != 360.)
2056 throw runtime_error(
"The range of *scat_aa_grid* must [0 360].");
2059 const Index stokes_dim = doit_scat_field.
ncols();
2060 assert(stokes_dim > 0 || stokes_dim < 5);
2070 if (atmosphere_dim == 1)
2073 (cloudbox_limits[1] - cloudbox_limits[0]) +1,
2074 1, 1, Nza, 1, stokes_dim));
2075 assert(
is_size(doit_scat_field,
2076 (cloudbox_limits[1] - cloudbox_limits[0]) +1,
2077 1, 1, scat_za_grid.
nelem(), 1, stokes_dim));
2079 else if (atmosphere_dim == 3)
2081 assert (
is_size( doit_i_field,
2082 (cloudbox_limits[1] - cloudbox_limits[0]) +1,
2083 (cloudbox_limits[3] - cloudbox_limits[2]) +1,
2084 (cloudbox_limits[5] - cloudbox_limits[4]) +1,
2085 Nza, Naa, stokes_dim));
2086 assert (
is_size( doit_scat_field,
2087 (cloudbox_limits[1] - cloudbox_limits[0]) +1,
2088 (cloudbox_limits[3] - cloudbox_limits[2]) +1,
2089 (cloudbox_limits[5] - cloudbox_limits[4]) +1,
2090 Nza, Naa, stokes_dim));
2095 os <<
"The atmospheric dimension must be 1D or 3D \n"
2096 <<
"for scattering calculations using the DOIT\n"
2097 <<
"module, but it is not. The value of *atmosphere_dim*\n"
2098 <<
"is " << atmosphere_dim <<
".";
2099 throw runtime_error( os.str() );
2102 if( !(doit_za_interp == 0 || doit_za_interp == 1 ) )
2103 throw runtime_error(
"Interpolation method is not defined. Use \n"
2104 "*doit_za_interpSet*.\n");
2106 if ( cloudbox_limits.
nelem()!= 2*atmosphere_dim)
2107 throw runtime_error(
2108 "*cloudbox_limits* is a vector which contains the"
2109 "upper and lower limit of the cloud for all "
2110 "atmospheric dimensions. So its dimension must"
2111 "be 2 x *atmosphere_dim*");
2113 if (doit_za_grid_size < 16)
2114 throw runtime_error(
2115 "*doit_za_grid_size* must be greater than 15 for"
2116 "accurate results");
2117 else if (doit_za_grid_size > 100)
2120 out1 <<
"Warning: doit_za_grid_size is very large which means that the \n"
2121 <<
"calculation will be very slow. The recommended value is 19.\n";
2127 Tensor4 pha_mat_local(doit_za_grid_size, scat_aa_grid.
nelem(),
2128 stokes_dim, stokes_dim, 0.);
2130 Tensor5 pha_mat_spt_local(pnd_field.
nbooks(), doit_za_grid_size,
2131 scat_aa_grid.
nelem(), stokes_dim, stokes_dim, 0.);
2135 nlinspace(za_grid, 0, 180, doit_za_grid_size);
2140 gridpos(gp_za_i, scat_za_grid, za_grid);
2142 Matrix itw_za_i(doit_za_grid_size, 2);
2146 Matrix doit_i_field_int(doit_za_grid_size, stokes_dim, 0);
2151 gridpos(gp_za, za_grid, scat_za_grid);
2157 Matrix doit_scat_field_org(doit_za_grid_size, stokes_dim, 0);
2162 grid_stepsize[0] = 180./(
Numeric)(doit_za_grid_size - 1);
2163 grid_stepsize[1] = 360./(
Numeric)(Naa - 1);
2165 Tensor3 product_field(doit_za_grid_size, Naa, stokes_dim, 0);
2167 if ( atmosphere_dim == 1 )
2169 Index scat_aa_index_local = 0;
2173 for (
Index p_index = 0;
2174 p_index <= cloudbox_limits[1]-cloudbox_limits[0];
2177 Numeric rtp_temperature_local =
2178 t_field(p_index + cloudbox_limits[0], 0, 0);
2180 for (
Index i = 0; i < stokes_dim; i++)
2182 if (doit_za_interp == 0)
2185 doit_i_field(p_index, 0, 0,
joker, 0, i), gp_za_i);
2187 else if (doit_za_interp == 1)
2190 for(
Index za = 0; za < za_grid.
nelem(); za++)
2192 doit_i_field_int(za, i) =
2194 doit_i_field(p_index, 0, 0,
joker, 0, i),
2204 for(
Index scat_za_index_local = 0;
2205 scat_za_index_local < doit_za_grid_size;
2206 scat_za_index_local++)
2209 Index index_zero = 0;
2212 out3 <<
"Calculate the phase matrix \n";
2214 scat_za_index_local,
2218 scat_aa_index_local,
2219 rtp_temperature_local,
2220 pha_mat_spt_agenda);
2223 pha_matCalc(pha_mat_local, pha_mat_spt_local, pnd_field,
2224 atmosphere_dim, p_index, 0,
2227 out3 <<
"Multiplication of phase matrix with incoming" <<
2234 for(
Index za_in = 0; za_in < doit_za_grid_size; za_in ++)
2236 for (
Index aa_in = 0; aa_in < Naa; ++ aa_in)
2241 for (
Index i = 0; i < stokes_dim; i++)
2243 for (
Index j = 0; j< stokes_dim; j++)
2245 product_field(za_in, aa_in, i) +=
2246 pha_mat_local(za_in, aa_in, i, j) *
2247 doit_i_field_int(za_in, j);
2254 out3 <<
"Compute integral. \n";
2255 for (
Index i = 0; i < stokes_dim; i++)
2257 doit_scat_field_org(scat_za_index_local, i)=
2267 for (
Index i = 0; i < stokes_dim; i++)
2269 if(doit_za_interp == 0)
2271 interp(doit_scat_field(p_index,
2278 doit_scat_field_org(
joker, i),
2283 for(
Index za = 0; za < scat_za_grid.
nelem(); za++)
2285 doit_scat_field(p_index, 0, 0, za, 0, i) =
2287 doit_scat_field_org(
joker, i),
2298 else if( atmosphere_dim == 3 ){
2300 for (
Index p_index = 0; p_index <= cloudbox_limits[1] - cloudbox_limits[0];
2303 for (
Index lat_index = 0; lat_index <=
2304 cloudbox_limits[3] - cloudbox_limits[2]; lat_index++)
2306 for (
Index lon_index = 0; lon_index <=
2307 cloudbox_limits[5] - cloudbox_limits[4]; lon_index++)
2310 Numeric rtp_temperature_local =
2311 t_field(p_index + cloudbox_limits[0],
2312 lat_index + cloudbox_limits[2],
2313 lon_index + cloudbox_limits[4]);
2316 for (
Index scat_aa_index_local = 1;
2317 scat_aa_index_local < Naa;
2318 scat_aa_index_local++)
2321 for (
Index i = 0; i < stokes_dim; i++)
2324 doit_i_field(p_index, lat_index, lon_index,
2325 joker, scat_aa_index_local, i), gp_za_i);
2328 for (
Index scat_za_index_local = 0;
2329 scat_za_index_local < doit_za_grid_size;
2330 scat_za_index_local++)
2333 out3 <<
"Calculate phase matrix \n";
2335 scat_za_index_local,
2339 scat_aa_index_local,
2340 rtp_temperature_local,
2341 pha_mat_spt_agenda);
2356 out3 <<
"Multiplication of phase matrix with" <<
2357 "incoming intensity \n";
2359 for(
Index za_in = 0; za_in < doit_za_grid_size; za_in ++)
2361 for (
Index aa_in = 0; aa_in < Naa; ++ aa_in)
2365 for (
Index i = 0; i < stokes_dim; i++)
2367 for (
Index j = 0; j< stokes_dim; j++)
2369 product_field(za_in, aa_in, i) +=
2370 pha_mat_local(za_in, aa_in, i, j) *
2371 doit_i_field_int(za_in, j);
2377 out3 <<
"Compute the integral \n";
2379 for (
Index i = 0; i < stokes_dim; i++)
2381 doit_scat_field_org(scat_za_index_local, i) =
2393 for (
Index i = 0; i < stokes_dim; i++)
2395 interp(doit_scat_field(p_index,
2399 scat_aa_index_local,
2402 doit_scat_field_org(
joker, i),
2412 out2 <<
" Finished scattered field.\n";
2418 Vector& doit_za_grid_opt,
2421 const Vector& scat_za_grid,
2422 const Index& doit_za_interp,
2433 chk_size(
"doit_i_field", doit_i_field,
2435 doit_i_field.
ncols());
2437 if(doit_i_field.
ncols()<1 || doit_i_field.
ncols()>4)
2438 throw runtime_error(
"The last dimension of *doit_i_field* corresponds\n"
2439 "to the Stokes dimension, therefore the number of\n"
2440 "columns in *doit_i_field* must be a number between\n"
2441 "1 and 4, but it is not!");
2443 if( !(doit_za_interp == 0 || doit_za_interp == 1 ) )
2444 throw runtime_error(
"Interpolation method is not defined. Use \n"
2445 "*doit_za_interpSet*.\n");
2447 if(scat_za_grid.
nelem() < 500)
2449 out1 <<
"Warning: The fine grid (*scat_za_grid*) has less than\n" <<
2450 "500 grid points which is likely not sufficient for\n" <<
2451 "grid_optimization\n" ;
2460 Matrix doit_i_field_opt_mat;
2461 doit_i_field_opt_mat = 0.;
2464 za_gridOpt(doit_za_grid_opt, doit_i_field_opt_mat,
2465 scat_za_grid, doit_i_field, acc,
2472 const Index& atmosphere_dim,
2479 if (atmosphere_dim != 1 && method ==
"polynomial")
2480 throw runtime_error(
2481 "Polynomial interpolation is only implemented for\n"
2482 "1D DOIT calculations as \n"
2483 "in 3D there can be numerical problems.\n"
2484 "Please use 'linear' interpolation method."
2487 if(method ==
"linear")
2489 else if (method ==
"polynomial")
2492 throw runtime_error(
"Possible interpolation methods are 'linear' "
2493 "and 'polynomial'.\n");
2504 Tensor4& doit_i_field1D_spectrum,
2505 const Index& atmfields_checked,
2506 const Index& atmgeom_checked,
2507 const Index& cloudbox_checked,
2508 const Index& cloudbox_on,
2510 const Agenda& doit_mono_agenda,
2511 const Index& doit_is_initialized,
2519 if( atmfields_checked != 1 )
2520 throw runtime_error(
"The atmospheric fields must be flagged to have "
2521 "passed a consistency check (atmfields_checked=1)." );
2522 if( atmgeom_checked != 1 )
2523 throw runtime_error(
"The atmospheric geometry must be flagged to have "
2524 "passed a consistency check (atmgeom_checked=1)." );
2525 if( cloudbox_checked != 1 )
2526 throw runtime_error(
"The cloudbox must be flagged to have "
2527 "passed a consistency check (cloudbox_checked=1)." );
2530 if (!cloudbox_on)
return;
2536 if( f_grid.
nelem() == 0 )
2537 throw runtime_error(
"The frequency grid is empty." );
2541 if (!doit_is_initialized)
2542 throw runtime_error(
2543 "Initialization method *DoitInit* has to be "
2545 "start of *ScatteringDoit*");
2553 Agenda l_doit_mono_agenda(doit_mono_agenda);
2562 for (
Index f_index = 0; f_index < nf; f_index ++)
2565 os <<
"Frequency: " << f_grid[f_index]/1e9 <<
" GHz \n" ;
2569 scat_i_lon, doit_i_field1D_spectrum,
2570 f_grid, f_index, l_doit_mono_agenda);
2580 Tensor4& doit_i_field1D_spectrum,
2584 const Index& f_index,
2588 const Vector& scat_za_grid,
2589 const Vector& scat_aa_grid,
2590 const Index& stokes_dim,
2591 const Index& atmosphere_dim,
2599 Index N_p = cloudbox_limits[1] - cloudbox_limits[0] +1;
2603 assert( f_index < f_grid.
nelem() );
2608 chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
2611 if (stokes_dim < 0 || stokes_dim > 4)
2612 throw runtime_error(
2613 "The dimension of stokes vector must be"
2616 if ( cloudbox_limits.nelem()!= 2*atmosphere_dim)
2617 throw runtime_error(
2618 "*cloudbox_limits* is a vector which contains the"
2619 "upper and lower limit of the cloud for all "
2620 "atmospheric dimensions. So its dimension must"
2621 "be 2 x *atmosphere_dim*");
2626 if(atmosphere_dim == 1)
2629 assert (
is_size( doit_i_field,
2630 (cloudbox_limits[1] - cloudbox_limits[0]) + 1,
2638 N_f, 2, 1, 1, N_za, 1, stokes_dim ));
2640 for (
Index za = 0; za < N_za; za++)
2642 for (
Index i = 0; i < stokes_dim; i++)
2646 scat_i_p(f_index, 0, 0, 0,
2648 doit_i_field(0, 0, 0, za, 0, i);
2650 scat_i_p(f_index, 1, 0, 0,
2652 doit_i_field(cloudbox_limits[1] - cloudbox_limits[0],
2656 doit_i_field1D_spectrum(f_index,
joker, za, i) =
2657 doit_i_field(
joker, 0, 0, za, 0, i);
2665 if(atmosphere_dim == 3)
2668 Index N_lat = cloudbox_limits[3] - cloudbox_limits[2] + 1;
2669 Index N_lon = cloudbox_limits[5] - cloudbox_limits[4] + 1;
2672 assert (
is_size( doit_i_field,
2673 cloudbox_limits[1] - cloudbox_limits[0] + 1,
2681 scat_i_p.
resize(N_f, 2, N_lat, N_lon, N_za, N_aa, stokes_dim);
2682 scat_i_lat.
resize(N_f, N_p, 2, N_lon, N_za, N_aa, stokes_dim);
2683 scat_i_lon.
resize(N_f, N_p, N_lat, 2, N_za, N_aa, stokes_dim);
2685 for (
Index za = 0; za < N_za; za++)
2687 for (
Index aa = 0; aa < N_aa; aa++)
2689 for (
Index i = 0; i < stokes_dim; i++)
2694 for (
Index lat = 0; lat < N_lat; lat++)
2696 for (
Index lon = 0; lon < N_lon; lon++)
2699 scat_i_p(f_index, 0, lat, lon,
2701 doit_i_field(0, lat, lon, za, aa, i);
2703 scat_i_p(f_index, 1, lat, lon,
2705 doit_i_field(cloudbox_limits[1]-cloudbox_limits[0],
2706 lat, lon, za, aa, i);
2712 for (
Index p = 0; p < N_p; p++)
2714 for (
Index lon = 0; lon < N_lon; lon++)
2717 scat_i_lat(f_index, p, 0, lon,
2719 doit_i_field(p, 0, lon, za, aa, i);
2721 scat_i_lat(f_index, p, 1, lon,
2723 doit_i_field(p, cloudbox_limits[3]-
2729 for (
Index lat = 0; lat < N_lat; lat++)
2732 scat_i_lon(f_index, p, lat, 0,
2734 doit_i_field(p, lat, 0, za, aa, i);
2736 scat_i_lon(f_index, p, lat, 1,
2738 doit_i_field(p, lat, cloudbox_limits[5]-
2739 cloudbox_limits[4], za, aa, i);
2756 const Index& atmfields_checked,
2757 const Index& atmgeom_checked,
2758 const Index& cloudbox_checked,
2759 const Index& doit_is_initialized,
2760 const Agenda& iy_main_agenda,
2761 const Index& atmosphere_dim,
2767 const Index& cloudbox_on,
2770 const Index& stokes_dim,
2772 const Agenda& blackbody_radiation_agenda,
2773 const Vector& scat_za_grid,
2774 const Vector& scat_aa_grid,
2775 const Index& rigorous,
2780 if( atmfields_checked != 1 )
2781 throw runtime_error(
"The atmospheric fields must be flagged to have "
2782 "passed a consistency check (atmfields_checked=1)." );
2783 if( atmgeom_checked != 1 )
2784 throw runtime_error(
"The atmospheric geometry must be flagged to have "
2785 "passed a consistency check (atmgeom_checked=1)." );
2786 if( cloudbox_checked != 1 )
2787 throw runtime_error(
"The cloudbox must be flagged to have "
2788 "passed a consistency check (cloudbox_checked=1)." );
2792 if (!cloudbox_on)
return;
2795 if (!doit_is_initialized)
2796 throw runtime_error(
2797 "Initialization method *DoitInit* has to be "
2798 "put before *CloudboxGetIncoming*");
2800 if( iy_unit !=
"1" ||
2804 os <<
"It is assumed that you use this method together with DOIT.\n"
2805 <<
"Usage of this method then demands that the *iy_main_agenda*\n"
2806 <<
"returns frequency based radiance (ie. [W/m2/Hz/sr]).\n"
2807 <<
"This requires that *iy_unit* is set to \"1\" and that\n"
2808 <<
"*blackbody_radiation_agenda uses *blackbody_radiationPlanck*\n"
2809 <<
"or a corresponding WSM.\n"
2810 <<
"At least one of these requirements is not met.\n";
2811 throw runtime_error( os.str() );
2815 Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
2817 Index Ni = stokes_dim;
2822 if( !(atmosphere_dim == 1 || atmosphere_dim == 3) )
2823 throw runtime_error(
"The atmospheric dimensionality must be 1 or 3.");
2824 if( scat_za_grid[0] != 0. || scat_za_grid[Nza-1] != 180. )
2825 throw runtime_error(
2826 "*scat_za_grid* must include 0 and 180 degrees as endpoints." );
2830 if( atmosphere_dim == 1 )
2833 scat_i_p.
resize( Nf, 2, 1, 1, Nza, 1, Ni );
2834 scat_i_lat.
resize( 0, 0, 0, 0, 0, 0, 0 );
2835 scat_i_lon.
resize( 0, 0, 0, 0, 0, 0, 0 );
2842 for (
Index boundary = 0; boundary <= 1; boundary++)
2844 pos[0] = z_field( cloudbox_limits[boundary], 0, 0 );
2848 los[0] = scat_za_grid[0];
2849 get_iy( ws, iy, t_field, z_field, vmr_field, 0, f_grid, pos, los,
2850 Vector(0), iy_main_agenda );
2851 scat_i_p(
joker, boundary, 0, 0, 0, 0,
joker ) = iy;
2853 for (
Index scat_za_index = 1; scat_za_index < Nza; scat_za_index ++)
2855 los[0] = scat_za_grid[scat_za_index];
2857 get_iy( ws, iy, t_field, z_field, vmr_field, 0, f_grid, pos, los,
2858 Vector(0), iy_main_agenda );
2860 scat_i_p(
joker, boundary, 0, 0, scat_za_index, 0,
joker ) = iy;
2864 for (
Index fi = 0; fi < Nf; fi ++)
2866 if( scat_i_p(fi,boundary,0,0,scat_za_index-1,0,0)/scat_i_p(fi,boundary,0,0,scat_za_index,0,0) > maxratio ||
2867 scat_i_p(fi,boundary,0,0,scat_za_index-1,0,0)/scat_i_p(fi,boundary,0,0,scat_za_index,0,0) < 1/maxratio )
2870 os <<
"ERROR: Radiance difference between interpolation\n"
2871 <<
"points is too large (factor " << maxratio <<
") to\n"
2872 <<
"safely interpolate. This might be due to za_grid\n"
2873 <<
"being too coarse or the radiance field being a\n"
2874 <<
"step-like function.\n";
2875 os <<
"Happens at boundary " << boundary <<
" between zenith\n"
2876 <<
"angels " << scat_za_grid[scat_za_index-1] <<
" and "
2877 << scat_za_grid[scat_za_index] <<
"deg for frequency"
2878 <<
"#" << fi <<
", where radiances are "
2879 << scat_i_p(fi,boundary,0,0,scat_za_index-1,0,0)
2880 <<
" and " << scat_i_p(fi,boundary,0,0,scat_za_index,0,0)
2881 <<
" W/(sr m2 Hz).";
2882 throw runtime_error(os.str());
2896 if( scat_aa_grid[0] != 0. || scat_aa_grid[Naa-1] != 360. )
2897 throw runtime_error(
2898 "*scat_aa_grid* must include 0 and 360 degrees as endpoints." );
2900 Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
2901 Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
2907 for(
Index i = 0; i<Naa; i++)
2908 aa_grid[i] = scat_aa_grid[i] - 180;
2911 scat_i_p.
resize( Nf, 2, Nlat_cloud, Nlon_cloud, Nza, Naa, Ni );
2912 scat_i_lat.
resize( Nf, Np_cloud, 2, Nlon_cloud, Nza, Naa, Ni );
2913 scat_i_lon.
resize( Nf, Np_cloud, Nlat_cloud, 2, Nza, Naa, Ni );
2921 for (
Index boundary = 0; boundary <= 1; boundary++)
2923 for (
Index lat_index = 0; lat_index < Nlat_cloud; lat_index++ )
2925 for (
Index lon_index = 0; lon_index < Nlon_cloud; lon_index++ )
2927 pos[2] = lon_grid[lon_index + cloudbox_limits[4]];
2928 pos[1] = lat_grid[lat_index + cloudbox_limits[2]];
2929 pos[0] = z_field(cloudbox_limits[boundary],
2930 lat_index + cloudbox_limits[2],
2931 lon_index + cloudbox_limits[4]);
2933 for (
Index scat_za_index = 0; scat_za_index < Nza;
2936 for (
Index scat_aa_index = 0; scat_aa_index < Naa;
2939 los[0] = scat_za_grid[scat_za_index];
2940 los[1] = aa_grid[scat_aa_index];
2945 if( ( scat_za_index != 0 &&
2946 scat_za_index != (Nza-1) ) ||
2947 scat_aa_index == 0 )
2949 get_iy( ws, iy, t_field, z_field, vmr_field, 0,
2950 f_grid, pos, los,
Vector(0),
2954 scat_i_p(
joker, boundary, lat_index, lon_index,
2955 scat_za_index, scat_aa_index,
joker) = iy;
2963 for (
Index boundary = 0; boundary <= 1; boundary++)
2965 for (
Index p_index = 0; p_index < Np_cloud; p_index++ )
2967 for (
Index lon_index = 0; lon_index < Nlon_cloud; lon_index++ )
2969 pos[2] = lon_grid[lon_index + cloudbox_limits[4]];
2970 pos[1] = lat_grid[cloudbox_limits[boundary+2]];
2971 pos[0] = z_field(p_index + cloudbox_limits[0],
2972 cloudbox_limits[boundary+2],
2973 lon_index + cloudbox_limits[4]);
2975 for (
Index scat_za_index = 0; scat_za_index < Nza;
2978 for (
Index scat_aa_index = 0; scat_aa_index < Naa;
2981 los[0] = scat_za_grid[scat_za_index];
2982 los[1] = aa_grid[scat_aa_index];
2986 if( ( scat_za_index != 0 &&
2987 scat_za_index != (Nza-1) ) ||
2988 scat_aa_index == 0 )
2990 get_iy( ws, iy, t_field, z_field, vmr_field, 0,
2991 f_grid, pos, los,
Vector(0),
2995 scat_i_lat(
joker, p_index, boundary, lon_index,
2996 scat_za_index, scat_aa_index,
joker) = iy;
3004 for (
Index boundary = 0; boundary <= 1; boundary++)
3006 for (
Index p_index = 0; p_index < Np_cloud; p_index++ )
3008 for (
Index lat_index = 0; lat_index < Nlat_cloud; lat_index++ )
3010 pos[2] = lon_grid[cloudbox_limits[boundary+4]];
3011 pos[1] = lat_grid[lat_index + cloudbox_limits[2]];
3012 pos[0] = z_field(p_index + cloudbox_limits[0],
3013 lat_index + cloudbox_limits[2],
3014 cloudbox_limits[boundary+4]);
3016 for (
Index scat_za_index = 0; scat_za_index < Nza;
3019 for (
Index scat_aa_index = 0; scat_aa_index < Naa;
3022 los[0] = scat_za_grid[scat_za_index];
3023 los[1] = aa_grid[scat_aa_index];
3027 if( ( scat_za_index != 0 &&
3028 scat_za_index != (Nza-1) ) ||
3029 scat_aa_index == 0 )
3031 get_iy( ws, iy, t_field, z_field, vmr_field, 0,
3032 f_grid, pos, los,
Vector(0),
3036 scat_i_lon(
joker, p_index, lat_index, boundary,
3037 scat_za_index, scat_aa_index,
joker) = iy;
3054 const Index& atmfields_checked,
3055 const Index& atmgeom_checked,
3056 const Index& cloudbox_checked,
3057 const Agenda& iy_main_agenda,
3058 const Index& atmosphere_dim,
3066 const Index& stokes_dim,
3068 const Agenda& blackbody_radiation_agenda,
3069 const Vector& scat_za_grid,
3070 const Vector& scat_aa_grid,
3074 if( atmfields_checked != 1 )
3075 throw runtime_error(
"The atmospheric fields must be flagged to have "
3076 "passed a consistency check (atmfields_checked=1)." );
3077 if( atmgeom_checked != 1 )
3078 throw runtime_error(
"The atmospheric geometry must be flagged to have "
3079 "passed a consistency check (atmgeom_checked=1)." );
3080 if( cloudbox_checked != 1 )
3081 throw runtime_error(
"The cloudbox must be flagged to have "
3082 "passed a consistency check (cloudbox_checked=1)." );
3085 if (!cloudbox_on)
return;
3088 if( iy_unit !=
"1" ||
3092 os <<
"It is assumed that you use this method together with DOIT.\n"
3093 <<
"Usage of this method then demands that the *iy_main_agenda*\n"
3094 <<
"returns frequency based radiance (ie. [W/m2/Hz/sr]).\n"
3095 <<
"This requires that *iy_unit* is set to \"1\" and that\n"
3096 <<
"*blackbody_radiation_agenda uses *blackbody_radiationPlanck*\n"
3097 <<
"or a corresponding WSM.\n"
3098 <<
"At least one of these requirements is not met.\n";
3099 throw runtime_error( os.str() );
3103 Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
3104 Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
3105 Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
3108 Index Ni = stokes_dim;
3113 if( atmosphere_dim != 3 )
3114 throw runtime_error(
"The atmospheric dimensionality must be 3.");
3115 if( scat_za_grid[0] != 0. || scat_za_grid[Nza-1] != 180. )
3116 throw runtime_error(
3117 "*scat_za_grid* must include 0 and 180 degrees as endpoints." );
3118 if( scat_aa_grid[0] != 0. || scat_aa_grid[Naa-1] != 360. )
3119 throw runtime_error(
3120 "*scat_aa_grid* must include 0 and 360 degrees as endpoints." );
3132 for(
Index i = 0; i<Naa; i++)
3133 aa_grid[i] = scat_aa_grid[i] - 180;
3140 scat_i_p.
resize(Nf, 2, Nlat_cloud, Nlon_cloud, Nza, Naa, Ni);
3141 scat_i_lat.
resize(Nf, Np_cloud, 2, Nlon_cloud, Nza, Naa, Ni);
3142 scat_i_lon.
resize(Nf, Np_cloud, Nlat_cloud, 2, Nza, Naa, Ni);
3148 pos[1] = lat_grid[cloudbox_limits[2]];
3149 pos[2] = lon_grid[cloudbox_limits[4]];
3150 los[1] = aa_grid[aa_index];
3153 for (
Index p_index = 0; p_index < Np_cloud; p_index++ )
3155 pos[0] = z_field( cloudbox_limits[0] + p_index, cloudbox_limits[2],
3156 cloudbox_limits[4] );
3158 for (
Index scat_za_index = 0; scat_za_index < Nza; scat_za_index++ )
3160 los[0] = scat_za_grid[scat_za_index];
3162 get_iy( ws, iy, t_field, z_field, vmr_field, 0, f_grid, pos, los,
3163 Vector(0), iy_main_agenda );
3165 for (
Index aa = 0; aa < Naa; aa ++)
3170 for (
Index lat = 0; lat < Nlat_cloud; lat ++)
3172 for (
Index lon = 0; lon < Nlon_cloud; lon ++)
3174 scat_i_p(
joker, 0, lat, lon, scat_za_index, aa,
3181 else if (p_index == Np_cloud-1)
3182 for (
Index lat = 0; lat < Nlat_cloud; lat ++)
3184 for (
Index lon = 0; lon < Nlon_cloud; lon ++)
3186 scat_i_p(
joker, 1, lat, lon, scat_za_index, aa,
3193 for (
Index lat = 0; lat < 2; lat ++)
3195 for (
Index lon = 0; lon < Nlon_cloud; lon ++)
3197 scat_i_lat(
joker, p_index, lat, lon,
3198 scat_za_index, aa,
joker)
3204 for (
Index lat = 0; lat < Nlat_cloud; lat ++)
3206 for (
Index lon = 0; lon < 2; lon ++)
3208 scat_i_lon(
joker, p_index, lat, lon,
3209 scat_za_index, aa,
joker)
3225 const Tensor4& doit_i_field1D_spectrum,
3228 const Index& jacobian_do,
3229 const Index& cloudbox_on,
3231 const Index& atmosphere_dim,
3236 const Index& stokes_dim,
3237 const Vector& scat_za_grid,
3238 const Vector& scat_aa_grid,
3240 const Index& rigorous,
3246 throw runtime_error(
3247 "This method does not provide any jacobians (jacobian_do must be 0)" );
3252 p_grid, lat_grid, lon_grid, z_field, rte_pos );
3256 doit_i_field1D_spectrum, gp_p, gp_lat, gp_lon,
3257 rte_los, cloudbox_on,
3258 cloudbox_limits, atmosphere_dim, stokes_dim,
3259 scat_za_grid, scat_aa_grid, f_grid,
"linear",
3270 const Tensor4& doit_i_field1D_spectrum,
3273 const Index& jacobian_do,
3274 const Index& cloudbox_on,
3276 const Index& atmosphere_dim,
3281 const Index& stokes_dim,
3282 const Vector& scat_za_grid,
3283 const Vector& scat_aa_grid,
3289 throw runtime_error(
3290 "This method does not provide any jacobians (jacobian_do must be 0)" );
3295 p_grid, lat_grid, lon_grid, z_field, rte_pos );
3299 doit_i_field1D_spectrum, gp_p, gp_lat, gp_lon,
3300 rte_los, cloudbox_on, cloudbox_limits,
3301 atmosphere_dim, stokes_dim,
3302 scat_za_grid, scat_aa_grid, f_grid,
"polynomial",
3311 const Tensor4& doit_i_field1D_spectrum,
3313 const Vector& scat_za_grid,
3315 const Index& f_index,
3316 const Index& atmosphere_dim,
3317 const Index& stokes_dim,
3322 if( atmosphere_dim!=1 )
3325 os <<
"This method is currently only implemented for 1D atmospheres!\n";
3326 throw runtime_error( os.str() );
3332 Index np = cloudbox_limits[1] - cloudbox_limits[0] +1;
3334 if( nf!=doit_i_field1D_spectrum.
nbooks() )
3337 os <<
"doit_i_field1D_spectrum has wrong size in frequency dimension.\n"
3338 << nf <<
" frequency points are expected, but doit_i_field1D_spectrum "
3339 <<
"contains " << doit_i_field1D_spectrum.
nbooks()
3340 <<
"frequency points.\n";
3341 throw runtime_error( os.str() );
3343 if( np!=doit_i_field1D_spectrum.
npages() )
3346 os <<
"doit_i_field1D_spectrum has wrong size in pressure level dimension.\n"
3347 << np <<
" pressure levels expected, but doit_i_field1D_spectrum "
3348 <<
"contains " << doit_i_field1D_spectrum.
npages()
3349 <<
"pressure levels.\n";
3350 throw runtime_error( os.str() );
3352 if( nza!=doit_i_field1D_spectrum.
nrows() )
3355 os <<
"doit_i_field1D_spectrum has wrong size in polar angle dimension.\n"
3356 << nza <<
" angles expected, but doit_i_field1D_spectrum "
3357 <<
"contains " << doit_i_field1D_spectrum.
nbooks()
3359 throw runtime_error( os.str() );
3361 if( stokes_dim!=doit_i_field1D_spectrum.
ncols() )
3364 os <<
"doit_i_field1D_spectrum has wrong stokes dimension.\n"
3365 <<
"Dimesnion " << stokes_dim
3366 <<
" expected, but doit_i_field1D_spectrum is dimesnion "
3367 << doit_i_field1D_spectrum.
nrows() <<
".\n";
3368 throw runtime_error( os.str() );
3372 doit_i_field.
resize(np,1,1,nza,1,stokes_dim);
3381 Index first_upwell = 0;
3382 assert(scat_za_grid[0]<90.);
3383 assert(scat_za_grid[scat_za_grid.
nelem()-1]>90.);
3384 while (scat_za_grid[first_upwell] < 90.)
3387 doit_i_field(np-1,0,0,
Range(0,first_upwell),0,
joker) =
3388 scat_i_p(f_index,1,0,0,
Range(0,first_upwell),0,
joker);
3390 if( cloudbox_limits[0] != 0 )
3391 doit_i_field(0,0,0,
Range(first_upwell,scat_za_grid.
nelem()-first_upwell),0,
joker) =
3392 scat_i_p(f_index,0,0,0,
Range(first_upwell,scat_za_grid.
nelem()-first_upwell),0,
joker);
3402 const Index& f_index,
3407 const Index& atmosphere_dim,
3409 const Index& all_frequencies,
3414 out2 <<
" Interpolate boundary clearsky field to obtain the initial field.\n";
3419 if(atmosphere_dim == 1)
3421 if(f_index == 0 || all_frequencies ==
true){
3423 if (f_grid.
nelem() != N_f){
3425 throw runtime_error(
" scat_i_p should have same frequency "
3426 " dimension as f_grid");
3430 throw runtime_error(
"scat_i_p should have exactly two elements "
3431 "in pressure dimension which correspond to the "
3432 "two cloudbox bounding pressure levels");
3442 doit_i_field.
resize((cloudbox_limits[1]- cloudbox_limits[0])+1,
3458 p_grid[
Range(cloudbox_limits[0],
3460 (cloudbox_limits[1]- cloudbox_limits[0]))],
3461 p_grid[
Range(cloudbox_limits[0],
3462 (cloudbox_limits[1]- cloudbox_limits[0])+1)]);
3464 Matrix itw((cloudbox_limits[1]- cloudbox_limits[0])+1, 2);
3469 for (
Index za_index = 0; za_index < N_za ; ++ za_index)
3471 for (
Index aa_index = 0; aa_index < N_aa ; ++ aa_index)
3473 for (
Index i = 0 ; i < N_i ; ++ i)
3511 else if(atmosphere_dim == 3)
3513 if (all_frequencies ==
false)
3514 throw runtime_error(
"Error in doit_i_fieldSetClearsky: For 3D "
3515 "all_frequencies option is not implemented \n");
3521 throw runtime_error(
" scat_i_p, scat_i_lat, scat_i_lon should have "
3522 "same frequency dimension");
3524 Index N_p = cloudbox_limits[1] - cloudbox_limits[0] + 1;
3527 throw runtime_error(
"scat_i_lat and scat_i_lon should have "
3528 "same pressure grid dimension as p_grid");
3531 Index N_lat = cloudbox_limits[3] - cloudbox_limits[2] + 1;
3533 if(scat_i_lon.
nshelves() != N_lat ||
3535 throw runtime_error(
"scat_i_p and scat_i_lon should have "
3536 "same latitude grid dimension as lat_grid");
3539 Index N_lon = cloudbox_limits[5] - cloudbox_limits[4] + 1;
3540 if(scat_i_lat.
nbooks() != N_lon ||
3541 scat_i_p.
nbooks() != N_lon ){
3542 throw runtime_error(
"scat_i_p and scat_i_lat should have "
3543 "same longitude grid dimension as lon_grid");
3546 throw runtime_error(
"scat_i_p should have exactly two elements "
3547 "in pressure dimension which correspond to the "
3548 "two cloudbox bounding pressure levels");
3552 throw runtime_error(
"scat_i_lat should have exactly two elements "
3553 "in latitude dimension which correspond to the "
3554 "two cloudbox bounding latitude levels");
3556 if(scat_i_lon.
nbooks() != 2){
3557 throw runtime_error(
"scat_i_lon should have exactly two elements "
3558 "in longitude dimension which correspond to the "
3559 "two cloudbox bounding longitude levels");
3562 if (scat_i_lat.
npages() != N_za ||
3563 scat_i_lon.
npages() != N_za){
3565 throw runtime_error(
" scat_i_p, scat_i_lat, scat_i_lon should have "
3566 "same zenith angle dimension");
3569 if (scat_i_lat.
nrows() != N_aa ||
3570 scat_i_lon.
nrows() != N_aa){
3572 throw runtime_error(
" scat_i_p, scat_i_lat, scat_i_lon should have "
3573 "same azimuth angle dimension");
3576 if (scat_i_lat.
ncols() != N_i ||
3577 scat_i_lon.
ncols() != N_i){
3579 throw runtime_error(
" scat_i_p, scat_i_lat, scat_i_lon should have "
3580 "same value for stokes_dim and can take only"
3581 "values 1,2,3 or 4");
3588 doit_i_field.
resize((cloudbox_limits[1]- cloudbox_limits[0])+1,
3589 (cloudbox_limits[3]- cloudbox_limits[2])+1,
3590 (cloudbox_limits[5]- cloudbox_limits[4])+1,
3598 ArrayOfGridPos lat_gp((cloudbox_limits[3]- cloudbox_limits[2])+1);
3599 ArrayOfGridPos lon_gp((cloudbox_limits[5]- cloudbox_limits[4])+1);
3606 p_grid[
Range(cloudbox_limits[0],
3608 (cloudbox_limits[1]- cloudbox_limits[0]))],
3609 p_grid[
Range(cloudbox_limits[0],
3610 (cloudbox_limits[1]- cloudbox_limits[0])+1)]);
3612 lat_grid[
Range(cloudbox_limits[2],
3614 (cloudbox_limits[3]- cloudbox_limits[2]))],
3615 lat_grid[
Range(cloudbox_limits[2],
3616 (cloudbox_limits[3]- cloudbox_limits[2])+1)]);
3618 lon_grid[
Range(cloudbox_limits[4],
3620 (cloudbox_limits[5]- cloudbox_limits[4]))],
3621 lon_grid[
Range(cloudbox_limits[4],
3622 (cloudbox_limits[5]- cloudbox_limits[4])+1)]);
3628 Matrix itw_p((cloudbox_limits[1]- cloudbox_limits[0])+1, 2);
3629 Matrix itw_lat((cloudbox_limits[3]- cloudbox_limits[2])+1, 2);
3630 Matrix itw_lon((cloudbox_limits[5]- cloudbox_limits[4])+1, 2);
3637 for (
Index lat_index = 0;
3638 lat_index <= (cloudbox_limits[3]-cloudbox_limits[2]); ++ lat_index)
3640 for (
Index lon_index = 0;
3641 lon_index <= (cloudbox_limits[5]-cloudbox_limits[4]);
3644 for (
Index za_index = 0; za_index < N_za ; ++ za_index)
3646 for (
Index aa_index = 0; aa_index < N_aa ; ++ aa_index)
3648 for (
Index i = 0 ; i < N_i ; ++ i)
3676 for (
Index p_index = 0;
3677 p_index <= (cloudbox_limits[1]-cloudbox_limits[0]) ; ++ p_index)
3679 for (
Index lon_index = 0;
3680 lon_index <= (cloudbox_limits[5]-cloudbox_limits[4]) ;
3683 for (
Index za_index = 0; za_index < N_za ; ++ za_index)
3685 for (
Index aa_index = 0; aa_index < N_aa ; ++ aa_index)
3687 for (
Index i = 0 ; i < N_i ; ++ i)
3692 za_index, aa_index, i);
3696 lon_index, za_index, aa_index, i);
3708 for (
Index p_index = 0;
3709 p_index <= (cloudbox_limits[1]-cloudbox_limits[0]); ++ p_index)
3711 for (
Index lat_index = 0;
3712 lat_index <= (cloudbox_limits[3]-cloudbox_limits[2]);
3715 for (
Index za_index = 0; za_index < N_za ; ++ za_index)
3717 for (
Index aa_index = 0; aa_index < N_aa ; ++ aa_index)
3719 for (
Index i = 0 ; i < N_i ; ++ i)
3722 VectorView target_field = doit_i_field(p_index,
3762 const Index& atmosphere_dim,
3763 const Index& stokes_dim,
3765 const Vector& doit_i_field_values,
3771 out2 <<
" Set initial field to constant values: " << doit_i_field_values <<
"\n";
3777 Index N_i = stokes_dim;
3778 Index Np_cloud = cloudbox_limits[1] - cloudbox_limits[0] + 1;
3779 Index Nlat_cloud = cloudbox_limits[3] - cloudbox_limits[2] + 1;
3780 Index Nlon_cloud = cloudbox_limits[5] - cloudbox_limits[4] + 1;
3785 chk_atm_grids( atmosphere_dim, p_grid, lat_grid, lon_grid );
3788 if (stokes_dim < 0 || stokes_dim > 4)
3789 throw runtime_error(
3790 "The dimension of stokes vector must be"
3793 if ( cloudbox_limits.nelem()!= 2*atmosphere_dim)
3794 throw runtime_error(
3795 "*cloudbox_limits* is a vector which contains the"
3796 "upper and lower limit of the cloud for all "
3797 "atmospheric dimensions. So its dimension must"
3798 "be 2 x *atmosphere_dim*");
3802 if(atmosphere_dim == 1)
3804 out3 <<
" atm_dim = 1\n";
3807 doit_i_field.
resize((cloudbox_limits[1] - cloudbox_limits[0])+1, 1, 1, N_za,
3812 for (
Index za_index = 0; za_index < N_za; za_index++)
3814 for (
Index i = 0; i < stokes_dim; i++)
3817 doit_i_field(cloudbox_limits[1]-cloudbox_limits[0], 0, 0, za_index,
3819 scat_i_p(0, 1, 0, 0, za_index, 0, i);
3821 doit_i_field(0, 0, 0, za_index, 0, i) =
3822 scat_i_p(0, 0, 0, 0, za_index, 0, i);
3823 for (
Index scat_p_index = 1; scat_p_index < cloudbox_limits[1] -
3824 cloudbox_limits[0]; scat_p_index++ )
3826 doit_i_field(scat_p_index, 0, 0, za_index, 0, i) = doit_i_field_values[i];
3832 if ( !
is_size(scat_i_p, 1, 2, Nlat_cloud,
3833 Nlon_cloud, N_za, N_aa, stokes_dim)
3834 || !
is_size(scat_i_lat, 1, Np_cloud, 2,
3835 Nlon_cloud, N_za, N_aa, stokes_dim)
3836 || !
is_size(scat_i_lon, 1, Np_cloud,
3837 Nlat_cloud, 2, N_za, N_aa, stokes_dim) )
3838 throw runtime_error(
3839 "One of the interface variables (*scat_i_p*, "
3840 "*scat_i_lat* or *scat_i_lon*) does not have "
3841 "the right dimensions. \n Probably you have "
3842 "calculated them before for another value of "
3848 out3 <<
"atm_dim = 3\n";
3849 doit_i_field.
resize((cloudbox_limits[1]- cloudbox_limits[0])+1,
3850 (cloudbox_limits[3]- cloudbox_limits[2])+1,
3851 (cloudbox_limits[5]- cloudbox_limits[4])+1,
3860 for (
Index za_index = 0; za_index < N_za; za_index++)
3862 for (
Index aa_index = 0; aa_index < N_aa; aa_index++)
3864 for (
Index i = 0; i < stokes_dim; i++)
3868 for (
Index lat_index = cloudbox_limits[2];
3869 lat_index <= cloudbox_limits[3]; lat_index++)
3871 for (
Index lon_index = cloudbox_limits[4];
3872 lon_index <= cloudbox_limits[5]; lon_index++)
3875 doit_i_field(cloudbox_limits[1]-cloudbox_limits[0],
3876 lat_index-cloudbox_limits[2],
3877 lon_index-cloudbox_limits[4],
3878 za_index, aa_index, i) =
3879 scat_i_p(0, 1, lat_index-cloudbox_limits[2],
3880 lon_index-cloudbox_limits[4],
3881 za_index, aa_index, i);
3883 doit_i_field(0, lat_index-cloudbox_limits[2],
3884 lon_index-cloudbox_limits[4],
3885 za_index, aa_index, i) =
3886 scat_i_p(0, 0, lat_index-cloudbox_limits[2],
3887 lon_index-cloudbox_limits[4],
3888 za_index, aa_index, i);
3892 for (
Index p_index = cloudbox_limits[0];
3893 p_index <= cloudbox_limits[1]; p_index++)
3897 for (
Index lon_index = cloudbox_limits[4];
3898 lon_index <= cloudbox_limits[5]; lon_index++)
3901 doit_i_field(p_index-cloudbox_limits[0],
3902 cloudbox_limits[3]-cloudbox_limits[2],
3903 lon_index-cloudbox_limits[4],
3904 za_index, aa_index, i) =
3905 scat_i_lat(0, p_index-cloudbox_limits[0],
3906 1, lon_index-cloudbox_limits[4],
3907 za_index, aa_index, i);
3909 doit_i_field(p_index-cloudbox_limits[0], 0,
3910 lon_index-cloudbox_limits[4],
3911 za_index, aa_index, i) =
3912 scat_i_lat(0, p_index-cloudbox_limits[0], 0,
3913 lon_index-cloudbox_limits[4],
3914 za_index, aa_index, i);
3918 for (
Index lat_index = cloudbox_limits[2];
3919 lat_index <= cloudbox_limits[3]; lat_index++)
3922 doit_i_field(p_index-cloudbox_limits[0],
3923 lat_index-cloudbox_limits[2],
3924 cloudbox_limits[5]-cloudbox_limits[4],
3925 za_index, aa_index, i) =
3926 scat_i_lon(0, p_index-cloudbox_limits[0],
3927 lat_index-cloudbox_limits[2], 1,
3928 za_index, aa_index, i);
3930 doit_i_field(p_index-cloudbox_limits[0],
3931 lat_index-cloudbox_limits[2],
3933 za_index, aa_index, i) =
3934 scat_i_lon(0, p_index-cloudbox_limits[0],
3935 lat_index-cloudbox_limits[2], 0,
3936 za_index, aa_index, i);
3943 for(
Index p_index = (cloudbox_limits[0]+1);
3944 p_index < cloudbox_limits[1] ;
3947 for (
Index lat_index = (cloudbox_limits[2]+1);
3948 lat_index < cloudbox_limits[3];
3951 for (
Index lon_index = (cloudbox_limits[4]+1);
3952 lon_index < cloudbox_limits[5];
3955 doit_i_field(p_index-cloudbox_limits[0],
3956 lat_index-cloudbox_limits[2],
3957 lon_index-cloudbox_limits[4],
3958 za_index, aa_index, i) =
3959 doit_i_field_values[i];