ARTS  2.2.66
ppath.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012 Patrick Eriksson <patrick.eriksson@chalmers.se>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA. */
17 
18 
19 
20 /*===========================================================================
21  === File description
22  ===========================================================================*/
23 
41 /*===========================================================================
42  === External declarations
43  ===========================================================================*/
44 
45 #include <cmath>
46 #include <stdexcept>
47 #include "agenda_class.h"
48 #include "array.h"
49 #include "arts_omp.h"
50 #include "auto_md.h"
51 #include "check_input.h"
52 #include "geodetic.h"
53 #include "math_funcs.h"
54 #include "messages.h"
55 #include "mystring.h"
56 #include "logic.h"
57 #include "poly_roots.h"
58 #include "ppath.h"
59 #include "refraction.h"
60 #include "rte.h"
61 #include "special_interp.h"
62 
63 extern const Numeric DEG2RAD;
64 extern const Numeric RAD2DEG;
65 
66 
67 
68 
69 
70 /*===========================================================================
71  === Precision variables
72  ===========================================================================*/
73 
74 // This variable defines the maximum allowed error tolerance for radius.
75 // The variable is, for example, used to check that a given a radius is
76 // consistent with the specified grid cell.
77 //
78 const Numeric RTOL = 1e-3;
79 
80 
81 // As RTOL but for latitudes and longitudes.
82 //
83 const Numeric LATLONTOL = 1e-8;
84 
85 
86 // Accuarcy for length comparisons.
87 //
88 const Numeric LACC = 1e-5;
89 
90 
91 // Values to apply if some calculation does not provide a solution
92 const Numeric R_NOT_FOUND = -1; // A value below zero
93 const Numeric L_NOT_FOUND = 99e99; // Some very large value for l/lat/lon
94 const Numeric LAT_NOT_FOUND = 99e99;
95 const Numeric LON_NOT_FOUND = 99e99;
96 
97 
98 
99 
100 
101 /*===========================================================================
102  === Functions related to geometrical propagation paths
103  ===========================================================================*/
104 
106 
118 Numeric geometrical_ppc( const Numeric& r, const Numeric& za )
119 {
120  assert( r > 0 );
121  assert( abs(za) <= 180 );
122 
123  return r * sin( DEG2RAD * abs(za) );
124 }
125 
126 
127 
129 
148  const Numeric& ppc,
149  const Numeric& a_za,
150  const Numeric& r )
151 {
152  assert( ppc >= 0 );
153  assert( abs(a_za) <= 180 );
154  assert( r >= ppc - RTOL );
155 
156  if( r > ppc )
157  {
158  Numeric za = RAD2DEG * asin( ppc / r );
159  if( abs(a_za) > 90 )
160  { za = 180 - za; }
161  if( a_za < 0 )
162  { za = -za; }
163  return za;
164  }
165  else
166  {
167  if( a_za > 0 )
168  { return 90; }
169  else
170  { return -90; }
171  }
172 }
173 
174 
175 
177 
191  const Numeric& ppc,
192  const Numeric& za )
193 {
194  assert( ppc >= 0 );
195  assert( abs(za) <= 180 );
196 
197  return ppc / sin( DEG2RAD * abs(za) );
198 }
199 
200 
201 
203 
219  const Numeric& za0,
220  const Numeric& lat0,
221  const Numeric& za )
222 {
223  assert( abs(za0) <= 180 );
224  assert( abs(za) <= 180 );
225  assert( ( za0 >= 0 && za >= 0 ) || ( za0 < 0 && za < 0 ) );
226 
227  return lat0 + za0 - za;
228 }
229 
230 
231 
233 
247  const Numeric& ppc,
248  const Numeric& r )
249 {
250  assert( ppc >= 0 );
251  assert( r >= ppc - RTOL );
252 
253  if( r > ppc )
254  { return sqrt( r*r - ppc*ppc ); }
255  else
256  { return 0; }
257 }
258 
259 
260 
262 
276  const Numeric& ppc,
277  const Numeric& l )
278 {
279  assert( ppc >= 0 );
280 
281  return sqrt( l*l + ppc*ppc );
282 }
283 
284 
285 
287 
300  const Numeric& ppc,
301  const Numeric& lat0,
302  const Numeric& za0,
303  const Numeric& lat )
304 {
305  assert( ppc >= 0 );
306  assert( abs(za0) <= 180 );
307  assert( ( za0 >= 0 && lat >= lat0 ) || ( za0 <= 0 && lat <= lat0 ) );
308 
309  // Zenith angle at the new latitude
310  const Numeric za = za0 + lat0 -lat;
311 
312  return geompath_r_at_za( ppc, za );
313 }
314 
315 
316 
318 
342  Vector& r,
343  Vector& lat,
344  Vector& za,
345  Numeric& lstep,
346  const Numeric& ppc,
347  const Numeric& r1,
348  const Numeric& lat1,
349  const Numeric& za1,
350  const Numeric& r2,
351  const bool& tanpoint,
352  const Numeric& lmax )
353 {
354  // Calculate length from tangent point, along the path for point 1 and 2.
355  // Length defined in the viewing direction, ie. negative at end of sensor
356  Numeric l1 = geompath_l_at_r( ppc, r1 );
357  if( abs(za1) > 90 )
358  { l1 *= -1; }
359  Numeric l2 = geompath_l_at_r( ppc, r2 );
360  if( l1 < 0 ) { l2 *= -1; }
361  if( tanpoint )
362  { l2 *= -1; }
363 
364  // Calculate needed number of steps, considering a possible length criterion
365  Index n;
366  if( lmax > 0 )
367  {
368  // The absolute value of the length distance is needed here
369  // We can't accept n=0, which is the case if l1 = l2
370  n = max( Index(1), Index( ceil( abs( l2 - l1 ) / lmax ) ) );
371  }
372  else
373  { n = 1; }
374 
375  // Length of path steps (note that lstep here can become negative)
376  lstep = ( l2 - l1 ) / (Numeric)n;
377 
378  // Allocate vectors and put in point 1
379  r.resize(n+1);
380  lat.resize(n+1);
381  za.resize(n+1);
382  r[0] = r1;
383  lat[0] = lat1;
384  za[0] = za1;
385 
386  // Loop steps (beside last) and calculate radius and zenith angle
387  for( Index i=1; i<n; i++ )
388  {
389  const Numeric l = l1 + lstep * (Numeric)i;
390  r[i] = geompath_r_at_l( ppc, l ); // Sign of l does not matter here
391  // Set a zenith angle to 80 or 100 depending on sign of l
392  za[i] = geompath_za_at_r( ppc, sign(za1)*(90-sign(l)*10), r[i] );
393  }
394 
395  // For maximum accuracy, set last radius to be exactly r2.
396  r[n] = r2; // Don't use r[n] below
397  za[n] = geompath_za_at_r( ppc, sign(za1)*(90-sign(l2)*10), r2 );
398 
399  // Ensure that zenith and nadir observations keep their zenith angle
400  if( abs(za1) < ANGTOL || abs(za1) > 180-ANGTOL )
401  { za = za1; }
402 
403  // Calculate latitudes
404  for( Index i=1; i<=n; i++ )
405  { lat[i] = geompath_lat_at_za( za1, lat1, za[i] ); }
406 
407  // Take absolute value of lstep
408  lstep = abs( lstep );
409 }
410 
411 
412 
413 
414 
415 /*===========================================================================
416  === Functions focusing on zenith and azimuth angles
417  ===========================================================================*/
418 
420 
444  Numeric& za,
445  Numeric& aa,
446  const Numeric& dx,
447  const Numeric& dy,
448  const Numeric& dz )
449 {
450  const Numeric r = sqrt( dx*dx + dy*dy + dz*dz );
451 
452  assert( r > 0 );
453 
454  za = RAD2DEG * acos( dz / r );
455  aa = RAD2DEG * atan2( dy, dx );
456 }
457 
458 
459 
461 
485  Numeric& dx,
486  Numeric& dy,
487  Numeric& dz,
488  const Numeric& za,
489  const Numeric& aa )
490 {
491  const Numeric zarad = DEG2RAD * za;
492  const Numeric aarad = DEG2RAD * aa;
493 
494  dz = cos( zarad );
495  dx = sin( zarad );
496  dy = sin( aarad ) * dx;
497  dx = cos( aarad ) * dx;
498 }
499 
500 
501 
503 
521  Matrix& R,
522  ConstVectorView vrot,
523  const Numeric& a )
524 {
525  assert( R.ncols() == 3 );
526  assert( R.nrows() == 3 );
527  assert( vrot.nelem() == 3 );
528 
529  const Numeric u = vrot[0];
530  const Numeric v = vrot[1];
531  const Numeric w = vrot[2];
532 
533  const Numeric u2 = u * u;
534  const Numeric v2 = v * v;
535  const Numeric w2 = w * w;
536 
537  assert( sqrt( u2 + v2 + w2 ) );
538 
539  const Numeric c = cos( DEG2RAD * a );
540  const Numeric s = sin( DEG2RAD * a );
541 
542  // Fill R
543  R(0,0) = u2 + (v2 + w2)*c;
544  R(0,1) = u*v*(1-c) - w*s;
545  R(0,2) = u*w*(1-c) + v*s;
546  R(1,0) = u*v*(1-c) + w*s;
547  R(1,1) = v2 + (u2+w2)*c;
548  R(1,2) = v*w*(1-c) - u*s;
549  R(2,0) = u*w*(1-c) - v*s;
550  R(2,1) = v*w*(1-c)+u*s;
551  R(2,2) = w2 + (u2+v2)*c;
552 }
553 
554 
555 
576 void map_daa(
577  Numeric& za,
578  Numeric& aa,
579  const Numeric& za0,
580  const Numeric& aa0,
581  const Numeric& aa_grid )
582 {
583  assert( abs( aa_grid ) <= 5 );
584 
585  Vector xyz(3);
586  Vector vrot(3);
587  Vector u(3);
588 
589  // Unit vector towards aa0 at za=90
590  //
591  zaaa2cart( xyz[0], xyz[1], xyz[2], 90, aa0 );
592 
593  // Find vector around which rotation shall be performed
594  //
595  // We can write this as cross([0 0 1],xyz). It turns out that the result
596  // of this operation is just [-y,x,0].
597  //
598  vrot[0] = -xyz[1];
599  vrot[1] = xyz[0];
600  vrot[2] = 0;
601 
602  // Unit vector towards aa0+aa at za=90
603  //
604  zaaa2cart( xyz[0], xyz[1], xyz[2], 90, aa0+aa_grid );
605 
606  // Apply rotation
607  //
608  Matrix R(3,3);
609  rotationmat3D( R, vrot, za0-90 );
610  mult( u, R, xyz );
611 
612  // Calculate za and aa for rotated u
613  //
614  cart2zaaa( za, aa, u[0], u[1], u[2] );
615 }
616 
617 
618 
619 
620 
621 /*===========================================================================
622  === Various functions
623  ===========================================================================*/
624 
626 
640  const Numeric& r,
641  const Numeric& za,
642  const Numeric& refr_index_air )
643 {
644  assert( r > 0 );
645  assert( abs(za) <= 180 );
646 
647  return r * refr_index_air * sin( DEG2RAD * abs(za) );
648 }
649 
650 
651 
653 
676  Numeric& lon,
677  const Numeric& lon5,
678  const Numeric& lon6 )
679 {
680  assert( lon6 >= lon5 );
681 
682  if( lon < lon5 && lon+180 <= lon6 )
683  { lon += 360; }
684  else if( lon > lon6 && lon-180 >= lon5 )
685  { lon -= 360; }
686 }
687 
688 
689 
691 
706  Index& it,
707  const Ppath ppath )
708 {
709  Numeric zmin = 99e99;
710  it = -1;
711  while( it < ppath.np-1 && ppath.pos(it+1,0) < zmin )
712  {
713  it++;
714  zmin = ppath.pos(it,0);
715  }
716  if( it == 0 || it == ppath.np-1 )
717  { it = -1; }
718 }
719 
720 
721 
722 
723 
724 /*===========================================================================
725  = 2D functions for surface and pressure level slope and tilt
726  ===========================================================================*/
727 
729 
744  const Numeric& lat1,
745  const Numeric& lat3,
746  const Numeric& r1,
747  const Numeric& r3,
748  const Numeric& lat )
749 {
750  return r1 + ( lat - lat1 ) * ( r3 - r1 ) / ( lat3 - lat1 );
751 }
752 
753 
754 
756 
783  Numeric& c1,
784  ConstVectorView lat_grid,
785  ConstVectorView refellipsoid,
786  ConstVectorView z_surf,
787  const GridPos& gp,
788  const Numeric& za )
789 {
790  Index i1 = gridpos2gridrange( gp, za >= 0 );
791  const Numeric r1 = refell2r( refellipsoid, lat_grid[i1] ) + z_surf[i1];
792  const Numeric r2 = refell2r( refellipsoid, lat_grid[i1+1] ) + z_surf[i1+1];
793  //
794  c1 = ( r2 - r1 ) / ( lat_grid[i1+1] - lat_grid[i1] );
795 }
796 
797 
798 
800 
818  Numeric& c1,
819  const Numeric& lat1,
820  const Numeric& lat2,
821  const Numeric& r1,
822  const Numeric& r2 )
823 {
824  c1 = ( r2 - r1 ) / ( lat2 - lat1 );
825 }
826 
827 
828 
829 
831 
845  const Numeric& r,
846  const Numeric& c1 )
847 {
848  // The tilt (in radians) is c1/r if c1 is converted to m/radian. So we get
849  // conversion RAD2DEG twice
850  return RAD2DEG * RAD2DEG * c1 / r;
851 }
852 
853 
854 
856 
876  const Numeric& za,
877  const Numeric& tilt )
878 {
879  assert( abs(za) <= 180 );
880 
881  // Yes, it shall be -tilt in both cases, if you wonder.
882  if( za > (90-tilt) || za < (-90-tilt) )
883  { return true; }
884  else
885  { return false; }
886 }
887 
888 
889 
891 
915  Numeric& lat,
916  Numeric& l,
917  const Numeric& r_hit,
918  const Numeric& r_start,
919  const Numeric& lat_start,
920  const Numeric& za_start,
921  const Numeric& ppc )
922 {
923  assert( abs( za_start ) <= 180 );
924  assert( r_start >= ppc );
925 
926  const Numeric absza = abs( za_start );
927 
928  // If above and looking upwards or r_hit below tangent point,
929  // we have no crossing:
930  if( ( r_start >= r_hit && absza <= 90 ) || ppc > r_hit )
931  { lat = LAT_NOT_FOUND; l = L_NOT_FOUND; }
932 
933  else
934  {
935  // Passages of tangent point
936  if( absza > 90 && r_start <= r_hit )
937  {
938  Numeric za = geompath_za_at_r( ppc, sign(za_start)*89, r_hit );
939  lat = geompath_lat_at_za( za_start, lat_start, za );
940  l = geompath_l_at_r( ppc, r_start ) + geompath_l_at_r( ppc, r_hit );
941  }
942 
943  else
944  {
945  Numeric za = geompath_za_at_r( ppc, za_start, r_hit );
946  lat = geompath_lat_at_za( za_start, lat_start, za );
947  l = abs( geompath_l_at_r( ppc, r_start ) -
948  geompath_l_at_r( ppc, r_hit ) );
949  assert( l > 0 );
950  }
951  }
952 }
953 
954 
955 
957 
993  const Numeric& rp,
994  const Numeric& za,
995  const Numeric& r0,
996  Numeric c1 )
997 {
998  const Numeric zaabs = abs(za);
999 
1000  assert( za != 0 );
1001  assert( zaabs != 180 );
1002  assert( abs(c1) > 0 ); // c1=0 should work, but unnecessary to use this func.
1003 
1004  // Convert slope to m/radian and consider viewing direction
1005  c1 *= RAD2DEG;
1006  if( za < 0 )
1007  { c1 = -c1; }
1008 
1009  // The nadir angle in radians, and cosine and sine of that angle
1010  const Numeric beta = DEG2RAD * ( 180 - zaabs );
1011  const Numeric cv = cos( beta );
1012  const Numeric sv = sin( beta );
1013 
1014  // Some repeated terms
1015  const Numeric r0s = r0*sv;
1016  const Numeric r0c = r0*cv;
1017  const Numeric cs = c1*sv;
1018  const Numeric cc = c1*cv;
1019 
1020  // The vector of polynomial coefficients
1021  //
1022  Index n = 6;
1023  Vector p0(n+1);
1024  //
1025  p0[0] = r0s - rp*sv;
1026  p0[1] = r0c + cs;
1027  p0[2] = -r0s/2 + cc;
1028  p0[3] = -r0c/6 - cs/2;
1029  p0[4] = r0s/24 - cc/6;
1030  p0[5] = r0c/120 + cs/24;
1031  p0[6] = -r0s/720 + cc/120;
1032  //
1033  // The accuracy when solving the polynomial equation gets worse when
1034  // approaching 0 and 180 deg. The solution is to let the start polynomial
1035  // order decrease when approaching these angles. The values below based on
1036  // practical experience, don't change without making extremely careful tests.
1037  //
1038  if( abs( 90 - zaabs ) > 89.9 )
1039  { n = 1; }
1040  else if( abs( 90 - zaabs ) > 75 )
1041  { n = 4; }
1042 
1043  // Calculate roots of the polynomial
1044  Matrix roots;
1045  int solutionfailure = 1;
1046  //
1047  while( solutionfailure)
1048  {
1049  roots.resize(n,2);
1050  Vector p;
1051  p = p0[Range(0,n+1)];
1052  solutionfailure = poly_root_solve( roots, p );
1053  if( solutionfailure )
1054  { n -= 1; assert( n > 0 ); }
1055  }
1056 
1057  // If r0=rp, numerical inaccuracy can give a false solution, very close
1058  // to 0, that we must throw away.
1059  Numeric dmin = 0;
1060  if( abs(r0-rp) < 1e-9 ) // 1 nm set based on practical experience.
1061  { dmin = 5e-12; }
1062 
1063  // Find the smallest root with imaginary part = 0, and real part > 0.
1064  Numeric dlat = 1.571; // Not interested in solutions above 90 deg!
1065  //
1066  for( Index i=0; i<n; i++ )
1067  {
1068  if( roots(i,1) == 0 && roots(i,0) > dmin && roots(i,0) < dlat )
1069  { dlat = roots(i,0); }
1070  }
1071 
1072  if( dlat < 1.57 ) // A somewhat smaller value than start one
1073  {
1074  // Convert back to degrees
1075  dlat = RAD2DEG * dlat;
1076 
1077  // Change sign if zenith angle is negative
1078  if( za < 0 )
1079  { dlat = -dlat; }
1080  }
1081  else
1082  { dlat = LAT_NOT_FOUND; }
1083 
1084  return dlat;
1085 }
1086 
1087 
1088 
1090 
1122  Numeric& r,
1123  Numeric& lat,
1124  Numeric& l,
1125  const Numeric& r_start0,
1126  const Numeric& lat_start,
1127  const Numeric& za_start,
1128  const Numeric& ppc,
1129  const Numeric& lat1,
1130  const Numeric& lat3,
1131  const Numeric& r1,
1132  const Numeric& r3,
1133  const bool& above )
1134 {
1135  const Numeric absza = abs( za_start );
1136 
1137  assert( absza <= 180 );
1138  assert( lat_start >=lat1 && lat_start <= lat3 );
1139 
1140  // Zenith case
1141  if( absza < ANGTOL )
1142  {
1143  if( above )
1144  { r = R_NOT_FOUND; lat = LAT_NOT_FOUND; l = L_NOT_FOUND; }
1145  else
1146  {
1147  lat = lat_start;
1148  r = rsurf_at_lat( lat1, lat3, r1, r3, lat );
1149  l = max( 1e-9, r - r_start0 ); // Max to ensure a small positive
1150  } // step, to handle numerical issues
1151  }
1152 
1153  // Nadir case
1154  else if( absza > 180-ANGTOL )
1155  {
1156  if( above )
1157  {
1158  lat = lat_start;
1159  r = rsurf_at_lat( lat1, lat3, r1, r3, lat );
1160  l = max( 1e-9, r_start0 - r ); // As above
1161  }
1162  else
1163  { r = R_NOT_FOUND; lat = LAT_NOT_FOUND; l = L_NOT_FOUND; }
1164  }
1165 
1166  // The general case
1167  else
1168  {
1169  const Numeric rmin = min( r1, r3 );
1170  const Numeric rmax = max( r1, r3 );
1171 
1172  // The case of negligible slope
1173  if( rmax-rmin < 1e-6 )
1174  {
1175  // Set r_start and r, considering impact of numerical problems
1176  Numeric r_start = r_start0;
1177  r = r1;
1178  if( above )
1179  { if( r_start < rmax ) { r_start = r = rmax; } }
1180  else
1181  { if( r_start > rmin ) { r_start = r = rmin; } }
1182 
1183  r_crossing_2d( lat, l, r, r_start, lat_start, za_start, ppc );
1184 
1185  // Check if inside [lat1,lat3]
1186  if( lat > lat3 || lat < lat1 )
1187  { r = R_NOT_FOUND; lat = LAT_NOT_FOUND; }
1188  }
1189 
1190  // With slope
1191  else
1192  {
1193  // Set r_start, considering impact of numerical problems
1194  Numeric r_start = r_start0;
1195  if( above )
1196  { if( r_start < rmin ) { r_start = rmin; } }
1197  else
1198  { if( r_start > rmax ) { r_start = rmax; } }
1199 
1200  Numeric za=999;
1201 
1202  // Calculate crossing with closest radius
1203  if( r_start > rmax )
1204  {
1205  r = rmax;
1206  r_crossing_2d( lat, l, r, r_start, lat_start, za_start, ppc );
1207  }
1208  else if( r_start < rmin )
1209  {
1210  r = rmin;
1211  r_crossing_2d( lat, l, r, r_start, lat_start, za_start, ppc );
1212  }
1213  else
1214  { r = r_start; lat = lat_start; l = 0; za = za_start; }
1215 
1216  // lat must be inside [lat1,lat3] if relevant to continue
1217  if( lat < lat1 || lat > lat3 )
1218  { r = R_NOT_FOUND; } // lat already set by r_crossing_2d
1219 
1220  // Otherwise continue from found point, considering the level slope
1221  else
1222  {
1223  // Level slope and radius at lat
1224  Numeric cpl;
1225  plevel_slope_2d( cpl, lat1, lat3, r1, r3 );
1226  const Numeric rpl = r1 + cpl*(lat-lat1);
1227 
1228  // Make adjustment if numerical problems
1229  if( above )
1230  { if( r < rpl ) { r = rpl; } }
1231  else
1232  { if( r > rpl ) { r = rpl; } }
1233 
1234  // Calculate zenith angle at (r,lat) (if <180, already set above)
1235  if( za > 180 ) // lat+za preserved (also with negative za)
1236  { za = lat_start + za_start - lat; };
1237 
1238  // Latitude distance from present point to actual crossing
1239  const Numeric dlat = rslope_crossing2d( r, za, rpl, cpl );
1240 
1241  // Update lat and check if still inside [lat1,lat3].
1242  // If yes, determine r and l
1243  lat += dlat;
1244  if( lat < lat1 || lat > lat3 )
1245  { r = R_NOT_FOUND; lat = LAT_NOT_FOUND; l = L_NOT_FOUND; }
1246  else
1247  {
1248  // It was tested to calculate r from geompath functions, but
1249  // appeared to give poorer accuracy around zenith/nadir
1250  r = rpl + cpl*dlat;
1251 
1252  // Zenith angle needed to check tangent point
1253  za = lat_start + za_start - lat;
1254 
1255  // Passage of tangent point requires special attention
1256  if( absza>90 && abs(za)<90 )
1257  { l = geompath_l_at_r( ppc, r_start ) +
1258  geompath_l_at_r( ppc, r ); }
1259  else
1260  { l = abs( geompath_l_at_r( ppc, r_start ) -
1261  geompath_l_at_r( ppc, r ) ); }
1262  }
1263  }
1264  }
1265  }
1266 }
1267 
1268 
1269 
1270 /*===========================================================================
1271  = 3D functions for level slope and tilt, and lat/lon crossings
1272  ===========================================================================*/
1273 
1275 
1295  const Numeric& lat1,
1296  const Numeric& lat3,
1297  const Numeric& lon5,
1298  const Numeric& lon6,
1299  const Numeric& r15,
1300  const Numeric& r35,
1301  const Numeric& r36,
1302  const Numeric& r16,
1303  const Numeric& lat,
1304  const Numeric& lon )
1305 {
1306  // We can't have any assert of *lat* and *lon* here as we can go outside
1307  // the ranges when called from *plevel_slope_3d*.
1308 
1309  if( lat == lat1 )
1310  { return r15 + ( lon - lon5 ) * ( r16 - r15 ) / ( lon6 -lon5 ); }
1311  else if( lat == lat3 )
1312  { return r35 + ( lon - lon5 ) * ( r36 - r35 ) / ( lon6 -lon5 ); }
1313  else if( lon == lon5 )
1314  { return r15 + ( lat - lat1 ) * ( r35 - r15 ) / ( lat3 -lat1 ); }
1315  else if( lon == lon6 )
1316  { return r16 + ( lat - lat1 ) * ( r36 - r16 ) / ( lat3 -lat1 ); }
1317  else
1318  {
1319  const Numeric fdlat = ( lat - lat1 ) / ( lat3 - lat1 );
1320  const Numeric fdlon = ( lon - lon5 ) / ( lon6 - lon5 );
1321  return (1-fdlat)*(1-fdlon)*r15 + fdlat*(1-fdlon)*r35 +
1322  (1-fdlat)*fdlon*r16 + fdlat*fdlon*r36;
1323  }
1324 }
1325 
1326 
1327 
1329 
1355  Numeric& c1,
1356  Numeric& c2,
1357  const Numeric& lat1,
1358  const Numeric& lat3,
1359  const Numeric& lon5,
1360  const Numeric& lon6,
1361  const Numeric& r15,
1362  const Numeric& r35,
1363  const Numeric& r36,
1364  const Numeric& r16,
1365  const Numeric& lat,
1366  const Numeric& lon,
1367  const Numeric& aa )
1368 {
1369  // Save time and avoid numerical problems if all r are equal
1370  if( r15==r35 && r15==r36 && r15==r16 && r35==r36 && r35==r16 && r36==r16 )
1371  {
1372  c1 = 0;
1373  c2 = 0;
1374  }
1375 
1376  else
1377  {
1378  // Radius at point of interest
1379  const Numeric r0 = rsurf_at_latlon( lat1, lat3, lon5, lon6,
1380  r15, r35, r36, r16, lat, lon );
1381 
1382  // Size of test angular distance. Unit is degrees.
1383  // dang = 1e-4 corresponds to about 10 m shift horisontally.
1384  // Smaller values should probably be avoided. For example, 1e-5 gave
1385  // significant c2 when it should be zero.
1386  const Numeric dang = 1e-4;
1387 
1388  Numeric lat2, lon2;
1389  latlon_at_aa( lat2, lon2, lat, lon, aa, dang );
1390  resolve_lon( lon2, lon5, lon6 );
1391  const Numeric dr1 = rsurf_at_latlon( lat1, lat3, lon5, lon6,
1392  r15, r35, r36, r16, lat2, lon2 ) - r0;
1393 
1394  latlon_at_aa( lat2, lon2, lat, lon, aa, 2*dang );
1395  resolve_lon( lon2, lon5, lon6 );
1396  const Numeric dr2 = rsurf_at_latlon( lat1, lat3, lon5, lon6,
1397  r15, r35, r36, r16, lat2, lon2 ) - r0;
1398 
1399  // Derive linear and quadratic coefficient
1400  c1 = 0.5 * ( 4*dr1 - dr2 );
1401  c2 = (dr1-c1)/(dang*dang);
1402  c1 /= dang;
1403  }
1404 }
1405 
1406 
1407 
1409 
1436  Numeric& c1,
1437  Numeric& c2,
1438  ConstVectorView lat_grid,
1439  ConstVectorView lon_grid,
1440  ConstVectorView refellipsoid,
1441  ConstMatrixView z_surf,
1442  const GridPos& gp_lat,
1443  const GridPos& gp_lon,
1444  const Numeric& aa )
1445 {
1446  Index ilat = gridpos2gridrange( gp_lat, abs( aa ) >= 0 );
1447  Index ilon = gridpos2gridrange( gp_lon, aa >= 0 );
1448 
1449  // Allow that we are at the upper edge of the grids only for special cases:
1450  const Index llat = lat_grid.nelem() - 1;
1451  // At North pole:
1452  if( ilat >= llat )
1453  {
1454  if( lat_grid[llat] > POLELAT )
1455  { ilat = llat - 1; }
1456  else
1457  { throw runtime_error( "The upper latitude end of the atmosphere "
1458  "reached, that is not allowed." ); }
1459  }
1460  // Complete 360deg lon coverage:
1461  if( ilon >= lon_grid.nelem() - 1 )
1462  {
1463  if( is_lon_cyclic( lon_grid ) )
1464  { ilon = 0; }
1465  else
1466  { throw runtime_error( "The upper longitude end of the atmosphere "
1467  "reached, that is not allowed." ); }
1468  }
1469 
1470  // Restore latitude and longitude values
1471  Vector itw(2);
1472  Numeric lat, lon;
1473  interpweights( itw, gp_lat );
1474  lat = interp( itw, lat_grid, gp_lat );
1475  interpweights( itw, gp_lon );
1476  lon = interp( itw, lon_grid, gp_lon );
1477 
1478  // Extract values that defines the grid cell
1479  const Numeric lat1 = lat_grid[ilat];
1480  const Numeric lat3 = lat_grid[ilat+1];
1481  const Numeric lon5 = lon_grid[ilon];
1482  const Numeric lon6 = lon_grid[ilon+1];
1483  const Numeric re1 = refell2r( refellipsoid, lat1 );
1484  const Numeric re3 = refell2r( refellipsoid, lat3 );
1485  const Numeric r15 = re1 + z_surf(ilat,ilon);
1486  const Numeric r35 = re3 + z_surf(ilat+1,ilon);
1487  const Numeric r36 = re3 + z_surf(ilat+1,ilon+1);
1488  const Numeric r16 = re1 + z_surf(ilat,ilon+1);
1489 
1490  plevel_slope_3d( c1, c2, lat1, lat3, lon5, lon6, r15, r35, r36, r16,
1491  lat, lon, aa );
1492 }
1493 
1494 
1495 
1497 
1512  const Numeric& rp,
1513  const Numeric& za,
1514  const Numeric& r0,
1515  Numeric c1,
1516  Numeric c2 )
1517 {
1518  // Even if the four corner radii of the grid cell differ, both c1 and c2 can
1519  // turn up to be 0. So in contrast to the 2D version, here we accept zero c.
1520 
1521  // Convert c1 and c2 from degrees to radians
1522  c1 *= RAD2DEG;
1523  c2 *= RAD2DEG*RAD2DEG;
1524 
1525  // The nadir angle in radians, and cosine and sine of that angle
1526  const Numeric beta = DEG2RAD * ( 180 - za );
1527  const Numeric cv = cos( beta );
1528  const Numeric sv = sin( beta );
1529 
1530  // Some repeated terms
1531  const Numeric r0s = r0*sv;
1532  const Numeric r0c = r0*cv;
1533  const Numeric c1s = c1*sv;
1534  const Numeric c1c = c1*cv;
1535  const Numeric c2s = c2*sv;
1536  const Numeric c2c = c2*cv;
1537 
1538  // The vector of polynomial coefficients
1539  //
1540  Index n = 6;
1541  Vector p0(n+1);
1542  //
1543  p0[0] = r0s - rp*sv;
1544  p0[1] = r0c + c1s;
1545  p0[2] = -r0s/2 + c1c + c2s;
1546  p0[3] = -r0c/6 - c1s/2 + c2c;
1547  p0[4] = r0s/24 - c1c/6 - c2s/2;
1548  p0[5] = r0c/120 + c1s/24 - c2c/6;
1549  p0[6] = -r0s/720 + c1c/120 + c2s/24;
1550  //
1551  // The accuracy when solving the polynomial equation gets worse when
1552  // approaching 0 and 180 deg. The solution is to let the start polynomial
1553  // order decrease when approaching these angles. The values below based
1554  // on practical experience, don't change without making extremly careful
1555  // tests.
1556  //
1557  if( abs( 90 - za ) > 89.9 )
1558  { n = 1; }
1559  else if( abs( 90 - za ) > 75 )
1560  { n = 4; }
1561 
1562  // Calculate roots of the polynomial
1563  Matrix roots;
1564  int solutionfailure = 1;
1565  //
1566  while( solutionfailure)
1567  {
1568  roots.resize(n,2);
1569  Vector p;
1570  p = p0[Range(0,n+1)];
1571  solutionfailure = poly_root_solve( roots, p );
1572  if( solutionfailure )
1573  { n -= 1; assert( n > 0 ); }
1574  }
1575 
1576  // If r0=rp, numerical inaccuracy can give a false solution, very close
1577  // to 0, that we must throw away.
1578  Numeric dmin = 0;
1579  if( r0 == rp )
1580  { dmin = 1e-6; }
1581 
1582  // Find the smallest root with imaginary part = 0, and real part > 0.
1583  //
1584  Numeric dlat = 1.571; // Not interested in solutions above 90 deg!
1585  //
1586  for( Index i=0; i<n; i++ )
1587  {
1588  if( roots(i,1) == 0 && roots(i,0) > dmin && roots(i,0) < dlat )
1589  { dlat = roots(i,0); }
1590  }
1591 
1592  if( dlat < 1.57 ) // A somewhat smaller value than start one
1593  {
1594  // Convert back to degrees
1595  dlat = RAD2DEG * dlat;
1596  }
1597  else
1598  { dlat = LAT_NOT_FOUND; }
1599 
1600  return dlat;
1601 }
1602 
1603 
1604 
1606 
1639  Numeric& lat,
1640  Numeric& lon,
1641  Numeric& l,
1642  const Numeric& r_hit,
1643  const Numeric& r_start,
1644  const Numeric& lat_start,
1645  const Numeric& lon_start,
1646  const Numeric& za_start,
1647  const Numeric& ppc,
1648  const Numeric& x,
1649  const Numeric& y,
1650  const Numeric& z,
1651  const Numeric& dx,
1652  const Numeric& dy,
1653  const Numeric& dz )
1654 {
1655  assert( za_start >= 0 );
1656  assert( za_start <= 180 );
1657 
1658  // If above and looking upwards or r_hit below tangent point,
1659  // we have no crossing:
1660  if( ( r_start >= r_hit && za_start <= 90 ) || ppc > r_hit )
1661  { lat = LAT_NOT_FOUND; lon = LON_NOT_FOUND; l = L_NOT_FOUND; }
1662 
1663  else
1664  {
1665  // Zenith/nadir
1666  if( za_start < ANGTOL || za_start > 180-ANGTOL )
1667  {
1668  l = abs( r_hit - r_start );
1669  lat = lat_start;
1670  lon = lon_start;
1671  }
1672 
1673  else
1674  {
1675  const Numeric p = x*dx + y*dy + z*dz;
1676  const Numeric pp = p * p;
1677  const Numeric q = x*x + y*y + z*z - r_hit*r_hit;
1678  const Numeric sq = sqrt( pp - q );
1679  const Numeric l1 = -p + sq;
1680  const Numeric l2 = -p - sq;
1681 
1682  const Numeric lmin = min( l1, l2 );
1683  const Numeric lmax = max( l1, l2 );
1684 
1685  // If r_start equals r_hit solutions just above zero can appear (that
1686  // shall be rejected). So we pick the biggest solution if lmin is
1687  // negative or just above zero.
1688  // (Tried to use "if( r_start != r_hit )", but failed occasionally)
1689  if( lmin < 1e-6 )
1690  { l = lmax; }
1691  else
1692  { l = lmin; }
1693  assert( l > 0 );
1694 
1695  lat = RAD2DEG * asin( ( z+dz*l ) / r_hit );
1696  lon = RAD2DEG * atan2( y+dy*l, x+dx*l );
1697  }
1698  }
1699 }
1700 
1701 
1702 
1703 
1704 
1705 /*===========================================================================
1706  === Basic functions for the Ppath structure
1707  ===========================================================================*/
1708 
1710 
1727  Ppath& ppath,
1728  const Index& atmosphere_dim,
1729  const Index& np )
1730 {
1731  assert( np > 0 );
1732  assert( atmosphere_dim >= 1 );
1733  assert( atmosphere_dim <= 3 );
1734 
1735  ppath.dim = atmosphere_dim;
1736  ppath.np = np;
1737  ppath.constant = -1;
1738 
1739  const Index npos = max( Index(2), atmosphere_dim );
1740  const Index nlos = max( Index(1), atmosphere_dim-1 );
1741 
1742  ppath.start_pos.resize( npos ); ppath.start_pos = -999;
1743  ppath.start_los.resize( nlos ); ppath.start_los = -999;
1744  ppath.start_lstep = 0;
1745  ppath.end_pos.resize( npos );
1746  ppath.end_los.resize( nlos );
1747  ppath.end_lstep = 0;
1748 
1749  ppath.pos.resize( np, npos );
1750  ppath.los.resize( np, nlos );
1751  ppath.r.resize( np );
1752  ppath.lstep.resize( np-1 );
1753 
1754  ppath.gp_p.resize( np );
1755  if( atmosphere_dim >= 2 )
1756  {
1757  ppath.gp_lat.resize( np );
1758  if( atmosphere_dim == 3 )
1759  { ppath.gp_lon.resize( np ); }
1760  }
1761 
1762  ppath_set_background( ppath, 0 );
1763  ppath.nreal.resize( np );
1764  ppath.ngroup.resize( np );
1765 }
1766 
1767 
1768 
1769 
1771 
1791  Ppath& ppath,
1792  const Index& case_nr )
1793 {
1794  switch ( case_nr )
1795  {
1796  case 0:
1797  ppath.background = "unvalid";
1798  break;
1799  case 1:
1800  ppath.background = "space";
1801  break;
1802  case 2:
1803  ppath.background = "surface";
1804  break;
1805  case 3:
1806  ppath.background = "cloud box level";
1807  break;
1808  case 4:
1809  ppath.background = "cloud box interior";
1810  break;
1811  case 9:
1812  ppath.background = "transmitter";
1813  break;
1814  default:
1815  ostringstream os;
1816  os << "Case number " << case_nr << " is not defined.";
1817  throw runtime_error(os.str());
1818  }
1819 }
1820 
1821 
1822 
1824 
1836 {
1837  if( ppath.background == "unvalid" )
1838  { return 0; }
1839  else if( ppath.background == "space" )
1840  { return 1; }
1841  else if( ppath.background == "surface" )
1842  { return 2; }
1843  else if( ppath.background == "cloud box level" )
1844  { return 3; }
1845  else if( ppath.background == "cloud box interior" )
1846  { return 4; }
1847  else if( ppath.background == "transmitter" )
1848  { return 9; }
1849  else
1850  {
1851  ostringstream os;
1852  os << "The string " << ppath.background
1853  << " is not a valid background case.";
1854  throw runtime_error(os.str());
1855  }
1856 }
1857 
1858 
1859 
1861 
1877  Ppath& ppath1,
1878  const Ppath& ppath2,
1879  const Index& ncopy )
1880 {
1881  Index n;
1882  if( ncopy < 0 )
1883  { n = ppath2.np; }
1884  else
1885  { n = ncopy; }
1886 
1887  assert( ppath1.np >= n );
1888 
1889  // The field np shall not be copied !
1890 
1891  ppath1.dim = ppath2.dim;
1892  ppath1.constant = ppath2.constant;
1893  ppath1.background = ppath2.background;
1894 
1895  // As index 0 always is included in the copying, the end point is covered
1896  ppath1.end_pos = ppath2.end_pos;
1897  ppath1.end_los = ppath2.end_los;
1898  ppath1.end_lstep = ppath2.end_lstep;
1899 
1900  // Copy start point only if copying fills up ppath1
1901  if( n == ppath1.np )
1902  {
1903  ppath1.start_pos = ppath2.start_pos;
1904  ppath1.start_los = ppath2.start_los;
1905  ppath1.start_lstep = ppath2.start_lstep;
1906  }
1907 
1908  ppath1.pos(Range(0,n),joker) = ppath2.pos(Range(0,n),joker);
1909  ppath1.los(Range(0,n),joker) = ppath2.los(Range(0,n),joker);
1910  ppath1.r[Range(0,n)] = ppath2.r[Range(0,n)];
1911  ppath1.nreal[Range(0,n)] = ppath2.nreal[Range(0,n)];
1912  ppath1.ngroup[Range(0,n)] = ppath2.ngroup[Range(0,n)];
1913  if( n > 1 )
1914  { ppath1.lstep[Range(0,n-1)] = ppath2.lstep[Range(0,n-1)]; }
1915 
1916  for( Index i=0; i<n; i++ )
1917  {
1918  gridpos_copy( ppath1.gp_p[i], ppath2.gp_p[i] );
1919 
1920  if( ppath1.dim >= 2 )
1921  { gridpos_copy( ppath1.gp_lat[i], ppath2.gp_lat[i] ); }
1922 
1923  if( ppath1.dim == 3 )
1924  { gridpos_copy( ppath1.gp_lon[i], ppath2.gp_lon[i] ); }
1925  }
1926 }
1927 
1928 
1929 
1931 
1949  Ppath& ppath1,
1950  const Ppath& ppath2 )
1951 {
1952  const Index n1 = ppath1.np;
1953  const Index n2 = ppath2.np;
1954 
1955  Ppath ppath;
1956  ppath_init_structure( ppath, ppath1.dim, n1 );
1957  ppath_copy( ppath, ppath1, -1 );
1958 
1959  ppath_init_structure( ppath1, ppath1.dim, n1 + n2 - 1 );
1960  ppath_copy( ppath1, ppath, -1 );
1961 
1962  // Append data from ppath2
1963  Index i1;
1964  for( Index i=1; i<n2; i++ )
1965  {
1966  i1 = n1 + i - 1;
1967 
1968  ppath1.pos(i1,0) = ppath2.pos(i,0);
1969  ppath1.pos(i1,1) = ppath2.pos(i,1);
1970  ppath1.los(i1,0) = ppath2.los(i,0);
1971  ppath1.r[i1] = ppath2.r[i];
1972  ppath1.nreal[i1] = ppath2.nreal[i];
1973  ppath1.ngroup[i1] = ppath2.ngroup[i];
1974  gridpos_copy( ppath1.gp_p[i1], ppath2.gp_p[i] );
1975 
1976  if( ppath1.dim >= 2 )
1977  { gridpos_copy( ppath1.gp_lat[i1], ppath2.gp_lat[i] ); }
1978 
1979  if( ppath1.dim == 3 )
1980  {
1981  ppath1.pos(i1,2) = ppath2.pos(i,2);
1982  ppath1.los(i1,1) = ppath2.los(i,1);
1983  gridpos_copy( ppath1.gp_lon[i1], ppath2.gp_lon[i] );
1984  }
1985 
1986  ppath1.lstep[i1-1] = ppath2.lstep[i-1];
1987  }
1988 
1989  if( ppath_what_background( ppath2 ) )
1990  { ppath1.background = ppath2.background; }
1991 
1992  ppath.start_pos = ppath2.start_pos;
1993  ppath.start_los = ppath2.start_los;
1994  ppath.start_lstep = ppath2.start_lstep;
1995 }
1996 
1997 
1998 
1999 
2000 
2001 /*===========================================================================
2002  === 1D/2D/3D start and end ppath functions
2003  ===========================================================================*/
2004 
2006 
2018  Numeric& r_start,
2019  Numeric& lat_start,
2020  Numeric& za_start,
2021  Index& ip,
2022  const Ppath& ppath )
2023 {
2024  // Number of points in the incoming ppath
2025  const Index imax = ppath.np - 1;
2026 
2027  // Extract starting radius, zenith angle and latitude
2028  r_start = ppath.r[imax];
2029  lat_start = ppath.pos(imax,1);
2030  za_start = ppath.los(imax,0);
2031 
2032  // Determine index of the pressure level being the lower limit for the
2033  // grid range of interest.
2034  //
2035  ip = gridpos2gridrange( ppath.gp_p[imax], za_start<=90 );
2036 }
2037 
2038 
2039 
2041 
2053  Ppath& ppath,
2054  ConstVectorView r_v,
2055  ConstVectorView lat_v,
2056  ConstVectorView za_v,
2057  ConstVectorView lstep,
2058  ConstVectorView n_v,
2059  ConstVectorView ng_v,
2060  ConstVectorView z_field,
2061  ConstVectorView refellipsoid,
2062  const Index& ip,
2063  const Index& endface,
2064  const Numeric& ppc )
2065 {
2066  // Number of path points
2067  const Index np = r_v.nelem();
2068 
2069  // Re-allocate ppath for return results and fill the structure
2070  ppath_init_structure( ppath, 1, np );
2071 
2072  ppath.constant = ppc;
2073 
2074  // Help variables that are common for all points.
2075  const Numeric r1 = refellipsoid[0] + z_field[ip];
2076  const Numeric dr = z_field[ip+1] - z_field[ip];
2077 
2078  for( Index i=0; i<np; i++ )
2079  {
2080  ppath.r[i] = r_v[i];
2081  ppath.pos(i,0) = r_v[i] - refellipsoid[0];
2082  ppath.pos(i,1) = lat_v[i];
2083  ppath.los(i,0) = za_v[i];
2084  ppath.nreal[i] = n_v[i];
2085  ppath.ngroup[i] = ng_v[i];
2086 
2087  ppath.gp_p[i].idx = ip;
2088  ppath.gp_p[i].fd[0] = ( r_v[i] - r1 ) / dr;
2089  ppath.gp_p[i].fd[1] = 1 - ppath.gp_p[i].fd[0];
2090  gridpos_check_fd( ppath.gp_p[i] );
2091 
2092  if( i > 0 )
2093  { ppath.lstep[i-1] = lstep[i-1]; }
2094  }
2095  gridpos_check_fd( ppath.gp_p[np-1] );
2096 
2097 
2098  //--- End point is the surface
2099  if( endface == 7 )
2100  { ppath_set_background( ppath, 2 ); }
2101 
2102  //--- End point is on top of a pressure level
2103  else if( endface <= 4 )
2104  { gridpos_force_end_fd( ppath.gp_p[np-1], z_field.nelem() ); }
2105 }
2106 
2107 
2108 
2110 
2122  Numeric& r_start,
2123  Numeric& lat_start,
2124  Numeric& za_start,
2125  Index& ip,
2126  Index& ilat,
2127  Numeric& lat1,
2128  Numeric& lat3,
2129  Numeric& r1a,
2130  Numeric& r3a,
2131  Numeric& r3b,
2132  Numeric& r1b,
2133  Numeric& rsurface1,
2134  Numeric& rsurface3,
2135  Ppath& ppath,
2136  ConstVectorView lat_grid,
2137  ConstMatrixView z_field,
2138  ConstVectorView refellipsoid,
2139  ConstVectorView z_surface
2140  )
2141 {
2142  // Number of points in the incoming ppath
2143  const Index imax = ppath.np - 1;
2144 
2145  // Extract starting radius, zenith angle and latitude
2146  r_start = ppath.r[imax];
2147  lat_start = ppath.pos(imax,1);
2148  za_start = ppath.los(imax,0);
2149 
2150  // Determine interesting latitude grid range and latitude end points of
2151  // the range.
2152  //
2153  ilat = gridpos2gridrange( ppath.gp_lat[imax], za_start >= 0 );
2154  //
2155  lat1 = lat_grid[ilat];
2156  lat3 = lat_grid[ilat+1];
2157 
2158  // Determine interesting pressure grid range. Do this first assuming that
2159  // the pressure levels are not tilted (that is, abs(za_start<=90) always
2160  // mean upward observation).
2161  // Set radius for the corners of the grid cell and the radial slope of
2162  // pressure level limits of the grid cell to match the found ip.
2163  //
2164  ip = gridpos2gridrange( ppath.gp_p[imax], abs(za_start) <= 90);
2165  //
2166  const Numeric re1 = refell2r( refellipsoid, lat_grid[ilat] );
2167  const Numeric re3 = refell2r( refellipsoid, lat_grid[ilat+1] );
2168  //
2169  r1a = re1 + z_field(ip,ilat); // lower-left
2170  r3a = re3 + z_field(ip,ilat+1); // lower-right
2171  r3b = re3 + z_field(ip+1,ilat+1); // upper-right
2172  r1b = re1 + z_field(ip+1,ilat); // upper-left
2173 
2174  // This part is a fix to catch start postions on top of a pressure level
2175  // that does not have an end fractional distance for the first step.
2176  {
2177  // Radius of lower and upper pressure level at the start position
2178  const Numeric rlow = rsurf_at_lat( lat1, lat3, r1a, r3a, lat_start );
2179  const Numeric rupp = rsurf_at_lat( lat1, lat3, r1b, r3b, lat_start );
2180  if( abs(r_start-rlow) < RTOL || abs(r_start-rupp) < RTOL )
2181  { gridpos_force_end_fd( ppath.gp_p[imax], z_field.nrows() ); }
2182  }
2183 
2184  // Slopes of pressure levels
2185  Numeric c2, c4;
2186  plevel_slope_2d( c2, lat1, lat3, r1a, r3a );
2187  plevel_slope_2d( c4, lat1, lat3, r1b, r3b );
2188 
2189  // Check if the LOS zenith angle happen to be between 90 and the zenith angle
2190  // of the pressure level (that is, 90 + tilt of pressure level), and in
2191  // that case if ip must be changed. This check is only needed when the
2192  // start point is on a pressure level.
2193  //
2194  if( is_gridpos_at_index_i( ppath.gp_p[imax], ip ) )
2195  {
2196  Numeric tilt = plevel_angletilt( r_start, c2 );
2197 
2198  if( is_los_downwards( za_start, tilt ) )
2199  {
2200  ip--;
2201  r1b = r1a; r3b = r3a;
2202  r1a = re1 + z_field(ip,ilat);
2203  r3a = re3 + z_field(ip,ilat+1);
2204  plevel_slope_2d( c2, lat1, lat3, r1a, r3a );
2205  }
2206  }
2207  else if( is_gridpos_at_index_i( ppath.gp_p[imax], ip+1 ) )
2208  {
2209  Numeric tilt = plevel_angletilt( r_start, c4 );
2210 
2211  if( !is_los_downwards( za_start, tilt ) )
2212  {
2213  ip++;
2214  r1a = r1b; r3a = r3b;
2215  r3b = re3 + z_field(ip+1,ilat+1);
2216  r1b = re1 + z_field(ip+1,ilat);
2217  plevel_slope_2d( c4, lat1, lat3, r1b, r3b );
2218  }
2219  }
2220 
2221  // Surface radius at latitude end points
2222  rsurface1 = re1 + z_surface[ilat];
2223  rsurface3 = re3 + z_surface[ilat+1];
2224 }
2225 
2226 
2227 
2229 
2241  Ppath& ppath,
2242  ConstVectorView r_v,
2243  ConstVectorView lat_v,
2244  ConstVectorView za_v,
2245  ConstVectorView lstep,
2246  ConstVectorView n_v,
2247  ConstVectorView ng_v,
2248  ConstVectorView lat_grid,
2249  ConstMatrixView z_field,
2250  ConstVectorView refellipsoid,
2251  const Index& ip,
2252  const Index& ilat,
2253  const Index& endface,
2254  const Numeric& ppc )
2255 {
2256  // Number of path points
2257  const Index np = r_v.nelem();
2258  const Index imax = np-1;
2259 
2260  // Re-allocate ppath for return results and fill the structure
2261  //
2262  ppath_init_structure( ppath, 2, np );
2263 
2264  ppath.constant = ppc;
2265 
2266  // Help variables that are common for all points.
2267  const Numeric dlat = lat_grid[ilat+1] - lat_grid[ilat];
2268  const Numeric z1low = z_field(ip,ilat);
2269  const Numeric z1upp = z_field(ip+1,ilat);
2270  const Numeric dzlow = z_field(ip,ilat+1) -z1low;
2271  const Numeric dzupp = z_field(ip+1,ilat+1) - z1upp;
2272  Numeric re = refell2r( refellipsoid, lat_grid[ilat] );
2273  const Numeric r1low = re + z1low;
2274  const Numeric r1upp = re + z1upp;
2275  re = refell2r( refellipsoid, lat_grid[ilat+1] );
2276  const Numeric drlow = re + z_field(ip,ilat+1) - r1low;
2277  const Numeric drupp = re + z_field(ip+1,ilat+1) - r1upp;
2278 
2279  for( Index i=0; i<np; i++ )
2280  {
2281  ppath.r[i] = r_v[i];
2282  ppath.pos(i,1) = lat_v[i];
2283  ppath.los(i,0) = za_v[i];
2284  ppath.nreal[i] = n_v[i];
2285  ppath.ngroup[i] = ng_v[i];
2286 
2287  // Weight in the latitude direction
2288  Numeric w = ( lat_v[i] - lat_grid[ilat] ) / dlat;
2289 
2290  // Radius of lower and upper face at present latitude
2291  const Numeric rlow = r1low + w * drlow;
2292  const Numeric rupp = r1upp + w * drupp;
2293 
2294  // Geometrical altitude of lower and upper face at present latitude
2295  const Numeric zlow = z1low + w * dzlow;
2296  const Numeric zupp = z1upp + w * dzupp;
2297 
2298  ppath.gp_p[i].idx = ip;
2299  ppath.gp_p[i].fd[0] = ( r_v[i] - rlow ) / ( rupp - rlow );
2300  ppath.gp_p[i].fd[1] = 1 - ppath.gp_p[i].fd[0];
2301  gridpos_check_fd( ppath.gp_p[i] );
2302 
2303  ppath.pos(i,0) = zlow + ppath.gp_p[i].fd[0] * ( zupp -zlow );
2304 
2305  ppath.gp_lat[i].idx = ilat;
2306  ppath.gp_lat[i].fd[0] = ( lat_v[i] - lat_grid[ilat] ) / dlat;
2307  ppath.gp_lat[i].fd[1] = 1 - ppath.gp_lat[i].fd[0];
2308  gridpos_check_fd( ppath.gp_lat[i] );
2309 
2310  if( i > 0 )
2311  { ppath.lstep[i-1] = lstep[i-1]; }
2312  }
2313  gridpos_check_fd( ppath.gp_p[imax] );
2314  gridpos_check_fd( ppath.gp_lat[imax] );
2315 
2316  // Do end-face specific tasks
2317  if( endface == 7 )
2318  { ppath_set_background( ppath, 2 ); }
2319 
2320  // Set fractional distance for end point
2321  //
2322  if( endface == 1 || endface == 3 )
2323  { gridpos_force_end_fd( ppath.gp_lat[np-1], lat_grid.nelem() ); }
2324  else if( endface == 2 || endface == 4 )
2325  { gridpos_force_end_fd( ppath.gp_p[np-1], z_field.nrows() ); }
2326 
2327  // Handle cases where exactly a corner is hit, or when slipping outside of
2328  // the grid box due to numerical inaccuarcy
2329  if( ppath.gp_p[imax].fd[0] < 0 || ppath.gp_p[imax].fd[1] < 0 )
2330  {
2331  gridpos_force_end_fd( ppath.gp_p[imax], z_field.nrows() );
2332  }
2333  if( ppath.gp_lat[imax].fd[0] < 0 || ppath.gp_lat[imax].fd[1] < 0 )
2334  {
2335  gridpos_force_end_fd( ppath.gp_lat[imax], lat_grid.nelem() );
2336  }
2337 }
2338 
2339 
2340 
2342 
2354  Numeric& r_start,
2355  Numeric& lat_start,
2356  Numeric& lon_start,
2357  Numeric& za_start,
2358  Numeric& aa_start,
2359  Index& ip,
2360  Index& ilat,
2361  Index& ilon,
2362  Numeric& lat1,
2363  Numeric& lat3,
2364  Numeric& lon5,
2365  Numeric& lon6,
2366  Numeric& r15a,
2367  Numeric& r35a,
2368  Numeric& r36a,
2369  Numeric& r16a,
2370  Numeric& r15b,
2371  Numeric& r35b,
2372  Numeric& r36b,
2373  Numeric& r16b,
2374  Numeric& rsurface15,
2375  Numeric& rsurface35,
2376  Numeric& rsurface36,
2377  Numeric& rsurface16,
2378  Ppath& ppath,
2379  ConstVectorView lat_grid,
2380  ConstVectorView lon_grid,
2381  ConstTensor3View z_field,
2382  ConstVectorView refellipsoid,
2383  ConstMatrixView z_surface )
2384 {
2385  // Index of last point in the incoming ppath
2386  const Index imax = ppath.np - 1;
2387 
2388  // Extract starting radius, zenith angle and latitude
2389  r_start = ppath.r[imax];
2390  lat_start = ppath.pos(imax,1);
2391  lon_start = ppath.pos(imax,2);
2392  za_start = ppath.los(imax,0);
2393  aa_start = ppath.los(imax,1);
2394 
2395  // Number of lat/lon
2396  const Index nlat = lat_grid.nelem();
2397  const Index nlon = lon_grid.nelem();
2398 
2399  // Lower index of lat and lon ranges of interest
2400  //
2401  // The longitude is undefined at the poles and as the azimuth angle
2402  // is defined in other way at the poles.
2403  //
2404  if( lat_start == 90 )
2405  {
2406  ilat = nlat - 2;
2407  GridPos gp_tmp;
2408  gridpos( gp_tmp, lon_grid, aa_start );
2409  if( aa_start < 180 )
2410  { ilon = gridpos2gridrange( gp_tmp, 1 ); }
2411  else
2412  { ilon = gridpos2gridrange( gp_tmp, 0 ); }
2413  }
2414  else if( lat_start == -90 )
2415  {
2416  ilat = 0;
2417  GridPos gp_tmp;
2418  gridpos( gp_tmp, lon_grid, aa_start );
2419  if( aa_start < 180 )
2420  { ilon = gridpos2gridrange( gp_tmp, 1 ); }
2421  else
2422  { ilon = gridpos2gridrange( gp_tmp, 0 ); }
2423  }
2424  else
2425  {
2426  if( lat_start > 0 )
2427  { ilat = gridpos2gridrange( ppath.gp_lat[imax], abs( aa_start )<90 ); }
2428  else
2429  { ilat = gridpos2gridrange( ppath.gp_lat[imax], abs( aa_start )<=90 );}
2430  if( lon_start < lon_grid[nlon-1] )
2431  { ilon = gridpos2gridrange( ppath.gp_lon[imax], aa_start >= 0 ); }
2432  else
2433  { ilon = nlon - 2; }
2434  }
2435  //
2436  lat1 = lat_grid[ilat];
2437  lat3 = lat_grid[ilat+1];
2438  lon5 = lon_grid[ilon];
2439  lon6 = lon_grid[ilon+1];
2440 
2441  // Determine interesting pressure grid range. Do this first assuming that
2442  // the pressure levels are not tilted (that is, abs(za_start<=90) always
2443  // mean upward observation).
2444  // Set radius for the corners of the grid cell and the radial slope of
2445  // pressure level limits of the grid cell to match the found ip.
2446  //
2447  ip = gridpos2gridrange( ppath.gp_p[imax], za_start <= 90 );
2448  //
2449  const Numeric re1 = refell2r( refellipsoid, lat_grid[ilat] );
2450  const Numeric re3 = refell2r( refellipsoid, lat_grid[ilat+1] );
2451  //
2452  r15a = re1 + z_field(ip,ilat,ilon);
2453  r35a = re3 + z_field(ip,ilat+1,ilon);
2454  r36a = re3 + z_field(ip,ilat+1,ilon+1);
2455  r16a = re1 + z_field(ip,ilat,ilon+1);
2456  r15b = re1 + z_field(ip+1,ilat,ilon);
2457  r35b = re3 + z_field(ip+1,ilat+1,ilon);
2458  r36b = re3 + z_field(ip+1,ilat+1,ilon+1);
2459  r16b = re1 + z_field(ip+1,ilat,ilon+1);
2460 
2461  // Check if the LOS zenith angle happen to be between 90 and the zenith angle
2462  // of the pressure level (that is, 90 + tilt of pressure level), and in
2463  // that case if ip must be changed. This check is only needed when the
2464  // start point is on a pressure level.
2465  //
2466  if( fabs(za_start-90) <= 10 ) // To save time. Ie. max tilt assumed =10 deg
2467  {
2468  if( is_gridpos_at_index_i( ppath.gp_p[imax], ip ) )
2469  {
2470  // Slope and angular tilt of lower pressure level
2471  Numeric c2a, c2b;
2472  plevel_slope_3d( c2a, c2b, lat1, lat3, lon5, lon6,
2473  r15a, r35a, r36a, r16a, lat_start, lon_start, aa_start );
2474  Numeric tilt = plevel_angletilt( r_start, c2a );
2475  // Negelect very small tilts, likely caused by numerical problems
2476  if( abs(tilt) > 1e-4 && is_los_downwards( za_start, tilt ) )
2477  {
2478  ip--;
2479  r15b = r15a; r35b = r35a; r36b = r36a; r16b = r16a;
2480  r15a = re1 + z_field(ip,ilat,ilon);
2481  r35a = re3 + z_field(ip,ilat+1,ilon);
2482  r36a = re3 + z_field(ip,ilat+1,ilon+1);
2483  r16a = re1 + z_field(ip,ilat,ilon+1);
2484  }
2485  }
2486  else if( is_gridpos_at_index_i( ppath.gp_p[imax], ip+1 ) )
2487  {
2488  // Slope and angular tilt of upper pressure level
2489  Numeric c4a, c4b;
2490  plevel_slope_3d( c4a, c4b, lat1, lat3 ,lon5, lon6,
2491  r15b, r35b, r36b, r16b, lat_start, lon_start, aa_start );
2492  Numeric tilt = plevel_angletilt( r_start, c4a );
2493  //
2494  if( !is_los_downwards( za_start, tilt ) )
2495  {
2496  ip++;
2497  r15a = r15b; r35a = r35b; r36a = r36b; r16a = r16b;
2498  r15b = re1 + z_field(ip+1,ilat,ilon);
2499  r35b = re3 + z_field(ip+1,ilat+1,ilon);
2500  r36b = re3 + z_field(ip+1,ilat+1,ilon+1);
2501  r16b = re1 + z_field(ip+1,ilat,ilon+1);
2502  }
2503  }
2504  }
2505 
2506  // Surface radius at latitude/longitude corner points
2507  rsurface15 = re1 + z_surface(ilat,ilon);
2508  rsurface35 = re3 + z_surface(ilat+1,ilon);
2509  rsurface36 = re3 + z_surface(ilat+1,ilon+1);
2510  rsurface16 = re1 + z_surface(ilat,ilon+1);
2511 }
2512 
2513 
2514 
2516 
2528  Ppath& ppath,
2529  ConstVectorView r_v,
2530  ConstVectorView lat_v,
2531  ConstVectorView lon_v,
2532  ConstVectorView za_v,
2533  ConstVectorView aa_v,
2534  ConstVectorView lstep,
2535  ConstVectorView n_v,
2536  ConstVectorView ng_v,
2537  ConstVectorView lat_grid,
2538  ConstVectorView lon_grid,
2539  ConstTensor3View z_field,
2540  ConstVectorView refellipsoid,
2541  const Index& ip,
2542  const Index& ilat,
2543  const Index& ilon,
2544  const Index& endface,
2545  const Numeric& ppc )
2546 {
2547  // Number of path points
2548  const Index np = r_v.nelem();
2549  const Index imax = np-1;
2550 
2551  // Re-allocate ppath for return results and fill the structure
2552  //
2553  ppath_init_structure( ppath, 3, np );
2554  //
2555  ppath.constant = ppc;
2556 
2557  // Help variables that are common for all points.
2558  const Numeric lat1 = lat_grid[ilat];
2559  const Numeric lat3 = lat_grid[ilat+1];
2560  const Numeric lon5 = lon_grid[ilon];
2561  const Numeric lon6 = lon_grid[ilon+1];
2562  const Numeric re1 = refell2r( refellipsoid, lat_grid[ilat] );
2563  const Numeric re3 = refell2r( refellipsoid, lat_grid[ilat+1] );
2564  const Numeric r15a = re1 + z_field(ip,ilat,ilon);
2565  const Numeric r35a = re3 + z_field(ip,ilat+1,ilon);
2566  const Numeric r36a = re3 + z_field(ip,ilat+1,ilon+1);
2567  const Numeric r16a = re1 + z_field(ip,ilat,ilon+1);
2568  const Numeric r15b = re1 + z_field(ip+1,ilat,ilon);
2569  const Numeric r35b = re3 + z_field(ip+1,ilat+1,ilon);
2570  const Numeric r36b = re3 + z_field(ip+1,ilat+1,ilon+1);
2571  const Numeric r16b = re1 + z_field(ip+1,ilat,ilon+1);
2572  const Numeric dlat = lat3 - lat1;
2573  const Numeric dlon = lon6 - lon5;
2574 
2575  for( Index i=0; i<np; i++ )
2576  {
2577  // Radius of pressure levels at present lat and lon
2578  const Numeric rlow = rsurf_at_latlon( lat1, lat3, lon5, lon6,
2579  r15a, r35a, r36a, r16a, lat_v[i], lon_v[i] );
2580  const Numeric rupp = rsurf_at_latlon( lat1, lat3, lon5, lon6,
2581  r15b, r35b, r36b, r16b, lat_v[i], lon_v[i] );
2582 
2583  // Position and LOS
2584  ppath.r[i] = r_v[i];
2585  ppath.pos(i,1) = lat_v[i];
2586  ppath.pos(i,2) = lon_v[i];
2587  ppath.los(i,0) = za_v[i];
2588  ppath.los(i,1) = aa_v[i];
2589  ppath.nreal[i] = n_v[i];
2590  ppath.ngroup[i] = ng_v[i];
2591 
2592  // Pressure grid index
2593  ppath.gp_p[i].idx = ip;
2594  ppath.gp_p[i].fd[0] = ( r_v[i] - rlow ) / ( rupp - rlow );
2595  ppath.gp_p[i].fd[1] = 1 - ppath.gp_p[i].fd[0];
2596  gridpos_check_fd( ppath.gp_p[i] );
2597 
2598  // Geometrical altitude
2599  const Numeric re = rsurf_at_latlon( lat1, lat3, lon5, lon6,
2600  re1, re3, re3, re1, lat_v[i], lon_v[i] );
2601  const Numeric zlow = rlow - re;
2602  const Numeric zupp = rupp - re;
2603  //
2604  ppath.pos(i,0) = zlow + ppath.gp_p[i].fd[0] * ( zupp -zlow );
2605 
2606  // Latitude grid index
2607  ppath.gp_lat[i].idx = ilat;
2608  ppath.gp_lat[i].fd[0] = ( lat_v[i] - lat1 ) / dlat;
2609  ppath.gp_lat[i].fd[1] = 1 - ppath.gp_lat[i].fd[0];
2610  gridpos_check_fd( ppath.gp_lat[i] );
2611 
2612  // Longitude grid index
2613  //
2614  // The longitude is undefined at the poles. The grid index is set to
2615  // the start point.
2616  //
2617  if( abs( lat_v[i] ) < POLELAT )
2618  {
2619  ppath.gp_lon[i].idx = ilon;
2620  ppath.gp_lon[i].fd[0] = ( lon_v[i] - lon5 ) / dlon;
2621  ppath.gp_lon[i].fd[1] = 1 - ppath.gp_lon[i].fd[0];
2622  gridpos_check_fd( ppath.gp_lon[i] );
2623  }
2624  else
2625  {
2626  ppath.gp_lon[i].idx = 0;
2627  ppath.gp_lon[i].fd[0] = 0;
2628  ppath.gp_lon[i].fd[1] = 1;
2629  }
2630 
2631  if( i > 0 )
2632  { ppath.lstep[i-1] = lstep[i-1]; }
2633  }
2634 
2635  // Do end-face specific tasks
2636  if( endface == 7 )
2637  { ppath_set_background( ppath, 2 ); }
2638 
2639  // Set fractional distance for end point
2640  //
2641  if( endface == 1 || endface == 3 )
2642  { gridpos_force_end_fd( ppath.gp_lat[imax], lat_grid.nelem() ); }
2643  else if( endface == 2 || endface == 4 )
2644  { gridpos_force_end_fd( ppath.gp_p[imax], z_field.npages() ); }
2645  else if( endface == 5 || endface == 6 )
2646  { gridpos_force_end_fd( ppath.gp_lon[imax], lon_grid.nelem() ); }
2647 
2648  // Handle cases where exactly a corner is hit, or when slipping outside of
2649  // the grid box due to numerical inaccuarcy
2650  if( ppath.gp_p[imax].fd[0] < 0 || ppath.gp_p[imax].fd[1] < 0 )
2651  {
2652  gridpos_force_end_fd( ppath.gp_p[imax], z_field.npages() );
2653  }
2654  if( ppath.gp_lat[imax].fd[0] < 0 || ppath.gp_lat[imax].fd[1] < 0 )
2655  {
2656  gridpos_force_end_fd( ppath.gp_lat[imax], lat_grid.nelem() );
2657  }
2658  if( ppath.gp_lon[imax].fd[0] < 0 || ppath.gp_lon[imax].fd[1] < 0 )
2659  {
2660  gridpos_force_end_fd( ppath.gp_lon[imax], lon_grid.nelem() );
2661  }
2662 }
2663 
2664 
2665 
2666 
2667 
2668 
2669 /*===========================================================================
2670  === Core functions for geometrical ppath_step calculations
2671  ===========================================================================*/
2672 
2674 
2700  Vector& r_v,
2701  Vector& lat_v,
2702  Vector& za_v,
2703  Numeric& lstep,
2704  Index& endface,
2705  const Numeric& r_start0,
2706  const Numeric& lat_start,
2707  const Numeric& za_start,
2708  const Numeric& ppc,
2709  const Numeric& lmax,
2710  const Numeric& ra,
2711  const Numeric& rb,
2712  const Numeric& rsurface )
2713 {
2714  Numeric r_start = r_start0;
2715 
2716  assert( rb > ra );
2717  assert( r_start >= ra - RTOL );
2718  assert( r_start <= rb + RTOL );
2719 
2720  // Shift radius if outside
2721  if( r_start < ra )
2722  { r_start = ra; }
2723  else if( r_start > rb )
2724  { r_start = rb; }
2725 
2726  Numeric r_end;
2727  bool tanpoint = false;
2728  endface = -1;
2729 
2730  // If upward, end point radius is always rb
2731  if( za_start <= 90 )
2732  { endface = 4; r_end = rb; }
2733 
2734  else
2735  {
2736  // Path reaches ra:
2737  if( ra > rsurface && ra > ppc )
2738  { endface = 2; r_end = ra; }
2739 
2740  // Path reaches the surface:
2741  else if( rsurface > ppc )
2742  { endface = 7; r_end = rsurface; }
2743 
2744  // The remaining option is a tangent point and back to rb
2745  else
2746  { endface = 4; r_end = rb; tanpoint = true; }
2747  }
2748 
2749  assert( endface > 0 );
2750 
2751  geompath_from_r1_to_r2( r_v, lat_v, za_v, lstep, ppc, r_start, lat_start,
2752  za_start, r_end, tanpoint, lmax );
2753 }
2754 
2755 
2756 
2758 
2784  Ppath& ppath,
2785  ConstVectorView z_field,
2786  ConstVectorView refellipsoid,
2787  const Numeric& z_surface,
2788  const Numeric& lmax )
2789 {
2790  // Starting radius, zenith angle and latitude
2791  Numeric r_start, lat_start, za_start;
2792 
2793  // Index of the pressure level being the lower limit for the
2794  // grid range of interest.
2795  Index ip;
2796 
2797  // Determine the variables defined above, and make asserts of input
2798  ppath_start_1d( r_start, lat_start, za_start, ip, ppath );
2799 
2800  // If the field "constant" is negative, this is the first call of the
2801  // function and the path constant shall be calculated.
2802  Numeric ppc;
2803  if( ppath.constant < 0 )
2804  { ppc = geometrical_ppc( r_start, za_start ); }
2805  else
2806  { ppc = ppath.constant; }
2807 
2808 
2809  // The path is determined by another function. Determine some variables
2810  // needed bý that function and call the function.
2811  //
2812  // Vars to hold found path points, path step length and coding for end face
2813  Vector r_v, lat_v, za_v;
2814  Numeric lstep;
2815  Index endface;
2816  //
2817  do_gridrange_1d( r_v, lat_v, za_v, lstep, endface,
2818  r_start, lat_start, za_start, ppc, lmax,
2819  refellipsoid[0]+z_field[ip], refellipsoid[0]+z_field[ip+1],
2820  refellipsoid[0]+z_surface );
2821 
2822  // Fill *ppath*
2823  const Index np = r_v.nelem();
2824  ppath_end_1d( ppath, r_v, lat_v, za_v, Vector(np-1,lstep), Vector(np,1),
2825  Vector(np,1), z_field, refellipsoid, ip, endface, ppc );
2826 }
2827 
2828 
2829 
2831 
2838  Vector& r_v,
2839  Vector& lat_v,
2840  Vector& za_v,
2841  Numeric& lstep,
2842  Index& endface,
2843  const Numeric& r_start0,
2844  const Numeric& lat_start0,
2845  const Numeric& za_start,
2846  const Numeric& l_start,
2847  const Index& icall,
2848  const Numeric& ppc,
2849  const Numeric& lmax,
2850  const Numeric& lat1,
2851  const Numeric& lat3,
2852  const Numeric& r1a,
2853  const Numeric& r3a,
2854  const Numeric& r3b,
2855  const Numeric& r1b,
2856  const Numeric& rsurface1,
2857  const Numeric& rsurface3 )
2858 {
2859  Numeric r_start = r_start0;
2860  Numeric lat_start = lat_start0;
2861 
2862  assert( icall < 4 );
2863 
2864  // Assert latitude and longitude
2865  assert( lat_start >= lat1 - LATLONTOL );
2866  assert( lat_start <= lat3 + LATLONTOL );
2867 
2868  // Shift latitude and longitude if outside
2869  if( lat_start < lat1 )
2870  { lat_start = lat1; }
2871  else if( lat_start > lat3 )
2872  { lat_start = lat3; }
2873 
2874  // Radius of lower and upper pressure level at the start position
2875  Numeric rlow = rsurf_at_lat( lat1, lat3, r1a, r3a, lat_start );
2876  Numeric rupp = rsurf_at_lat( lat1, lat3, r1b, r3b, lat_start );
2877 
2878  // Assert radius (some extra tolerance is needed for radius)
2879  assert( r_start >= rlow - RTOL );
2880  assert( r_start <= rupp + RTOL );
2881 
2882  // Shift radius if outside
2883  if( r_start < rlow )
2884  { r_start = rlow; }
2885  else if( r_start > rupp )
2886  { r_start = rupp; }
2887 
2888  // Position and LOS in cartesian coordinates
2889  Numeric x, z, dx, dz;
2890  poslos2cart( x, z, dx, dz, r_start, lat_start, za_start );
2891 
2892  // Some booleans to check for recursive call
2893  bool unsafe = false;
2894  bool do_surface = false;
2895 
2896  // Determine the position of the end point
2897  //
2898  endface = 0;
2899  //
2900  Numeric r_end, lat_end, l_end;
2901 
2902  // Zenith and nadir looking are handled as special cases
2903  const Numeric absza = abs( za_start );
2904 
2905  // Zenith looking
2906  if( absza < ANGTOL )
2907  {
2908  l_end = rupp - r_start;
2909  endface = 4;
2910  }
2911 
2912  // Nadir looking
2913  else if( absza > 180-ANGTOL )
2914  {
2915  const Numeric rsurface = rsurf_at_lat( lat1, lat3, rsurface1, rsurface3,
2916  lat_start );
2917 
2918  if( rlow > rsurface )
2919  {
2920  l_end = r_start - rlow;
2921  endface = 2;
2922  }
2923  else
2924  {
2925  l_end = r_start - rsurface;
2926  endface = 7;
2927  }
2928  }
2929 
2930  else
2931  {
2932  unsafe = true;
2933 
2934  // Calculate correction terms for the position to compensate for
2935  // numerical inaccuracy.
2936  //
2937  Numeric r_corr, lat_corr;
2938  //
2939  cart2pol( r_corr, lat_corr, x, z, lat_start, za_start );
2940  //
2941  r_corr -= r_start;
2942  lat_corr -= lat_start;
2943 
2944  // The end point is found by testing different step lengths until the
2945  // step length has been determined by a precision of *l_acc*.
2946  //
2947  // The first step is to found a length that goes outside the grid cell,
2948  // to find an upper limit. The lower limit is set to 0. The upper
2949  // limit is either set to l_start or it is estimated from the size of
2950  // the grid cell.
2951  // The search algorith is bisection, the next length to test is the
2952  // mean value of the minimum and maximum length limits.
2953  //
2954  if( l_start > 0 )
2955  { l_end = l_start; }
2956  else
2957  { l_end = 2 * ( rupp - rlow ); }
2958  //
2959  Numeric l_in = 0, l_out = l_end;
2960  bool ready = false, startup = true;
2961 
2962  if( rsurface1+RTOL >= r1a || rsurface3+RTOL >= r3a )
2963  { do_surface = true; }
2964 
2965  while( !ready )
2966  {
2967  cart2pol( r_end, lat_end, x+dx*l_end, z+dz*l_end, lat_start,
2968  za_start );
2969  r_end -= r_corr;
2970  lat_end -= lat_corr;
2971 
2972  bool inside = true;
2973 
2974  rlow = rsurf_at_lat( lat1, lat3, r1a, r3a, lat_end );
2975 
2976  if( do_surface )
2977  {
2978  const Numeric r_surface = rsurf_at_lat( lat1, lat3, rsurface1,
2979  rsurface3, lat_end );
2980  if( r_surface >= rlow && r_end <= r_surface )
2981  { inside = false; endface = 7; }
2982  }
2983 
2984  if( inside )
2985  {
2986  if( lat_end < lat1 )
2987  { inside = false; endface = 1; }
2988  else if( lat_end > lat3 )
2989  { inside = false; endface = 3; }
2990  else if( r_end < rlow )
2991  { inside = false; endface = 2; }
2992  else
2993  {
2994  rupp = rsurf_at_lat( lat1, lat3, r1b, r3b, lat_end );
2995  if( r_end > rupp )
2996  { inside = false; endface = 4; }
2997  }
2998  }
2999 
3000  if( startup )
3001  {
3002  if( inside )
3003  {
3004  l_in = l_end;
3005  l_end *= 5;
3006  }
3007  else
3008  {
3009  l_out = l_end;
3010  l_end = ( l_out + l_in ) / 2;
3011  startup = false;
3012  }
3013  }
3014  else
3015  {
3016  if( inside )
3017  { l_in = l_end; }
3018  else
3019  { l_out = l_end; }
3020 
3021  if( ( l_out - l_in ) < LACC )
3022  { ready = true; }
3023  else
3024  { l_end = ( l_out + l_in ) / 2; }
3025  }
3026  }
3027  }
3028 
3029  //--- Create return vectors
3030  //
3031  Index n = 1;
3032  //
3033  if( lmax > 0 )
3034  {
3035  n = Index( ceil( abs( l_end / lmax ) ) );
3036  if( n < 1 )
3037  { n = 1; }
3038  }
3039  //
3040  r_v.resize( n+1 );
3041  lat_v.resize( n+1 );
3042  za_v.resize( n+1 );
3043  //
3044  r_v[0] = r_start;
3045  lat_v[0] = lat_start;
3046  za_v[0] = za_start;
3047  //
3048  lstep = l_end / (Numeric)n;
3049  Numeric l;
3050  bool ready = true;
3051  //
3052  for( Index j=1; j<=n; j++ )
3053  {
3054  l = lstep * (Numeric)j;
3055  cart2poslos( r_v[j], lat_v[j], za_v[j], x+dx*l, z+dz*l, dx, dz, ppc,
3056  lat_start, za_start );
3057 
3058  if( j < n )
3059  {
3060  if( unsafe )
3061  {
3062  // Check that r_v[j] is above lower pressure level and the
3063  // surface. This can fail around tangent points. For p-levels
3064  // with constant r this is easy to handle analytically, but the
3065  // problem is tricky in the general case with a non-spherical
3066  // geometry, and this crude solution is used instead. Not the
3067  // most elegant solution, but it works! Added later the same
3068  // check for upper level, after getting assert in that direction.
3069  // The z_field was crazy, but still formerly correct.
3070  rlow = rsurf_at_lat( lat1, lat3, r1a, r3a, lat_v[j] );
3071  if( do_surface )
3072  {
3073  const Numeric r_surface = rsurf_at_lat( lat1, lat3,
3074  rsurface1, rsurface3, lat_v[j] );
3075  const Numeric r_test = max( r_surface, rlow );
3076  if( r_v[j] < r_test )
3077  { ready = false; break; }
3078  }
3079  else if( r_v[j] < rlow )
3080  { ready = false; break; }
3081 
3082  rupp = rsurf_at_lat( lat1, lat3, r1b, r3b, lat_v[j] );
3083  if( r_v[j] > rupp )
3084  { ready = false; break; }
3085  }
3086  }
3087  else // j==n
3088  {
3089  if( unsafe )
3090  {
3091  // Set end point to be consistent with found endface.
3092  //
3093  if( endface == 1 )
3094  { lat_v[n] = lat1; }
3095  else if( endface == 2 )
3096  { r_v[n] = rsurf_at_lat( lat1, lat3, r1a, r3a, lat_v[n] ); }
3097  else if( endface == 3 )
3098  { lat_v[n] = lat3; }
3099  else if( endface == 4 )
3100  { r_v[n] = rsurf_at_lat( lat1, lat3, r1b, r3b, lat_v[n] ); }
3101  else if( endface == 7 )
3102  { r_v[n] = rsurf_at_lat( lat1, lat3, rsurface1, rsurface3,
3103  lat_v[n] ); }
3104  }
3105  }
3106  }
3107 
3108  if( !ready )
3109  { // If an "outside" point found, restart with l as start search length
3110  do_gridcell_2d_byltest( r_v, lat_v, za_v, lstep, endface, r_start,
3111  lat_start, za_start, l, icall+1, ppc, lmax,
3112  lat1, lat3, r1a, r3a, r3b, r1b,
3113  rsurface1, rsurface3 );
3114  }
3115 }
3116 
3117 
3118 
3120 
3137  Ppath& ppath,
3138  ConstVectorView lat_grid,
3139  ConstMatrixView z_field,
3140  ConstVectorView refellipsoid,
3141  ConstVectorView z_surface,
3142  const Numeric& lmax )
3143 {
3144  // Radius, zenith angle and latitude of start point.
3145  Numeric r_start, lat_start, za_start;
3146 
3147  // Lower grid index for the grid cell of interest.
3148  Index ip, ilat;
3149 
3150  // Radii and latitudes set by *ppath_start_2d*.
3151  Numeric lat1, lat3, r1a, r3a, r3b, r1b, rsurface1, rsurface3;
3152 
3153  // Determine the variables defined above and make all possible asserts
3154  ppath_start_2d( r_start, lat_start, za_start, ip, ilat,
3155  lat1, lat3, r1a, r3a, r3b, r1b, rsurface1, rsurface3,
3156  ppath, lat_grid, z_field, refellipsoid, z_surface );
3157 
3158  // If the field "constant" is negative, this is the first call of the
3159  // function and the path constant shall be calculated.
3160  Numeric ppc;
3161  if( ppath.constant < 0 )
3162  { ppc = geometrical_ppc( r_start, za_start ); }
3163  else
3164  { ppc = ppath.constant; }
3165 
3166  // Vars to hold found path points, path step length and coding for end face
3167  Vector r_v, lat_v, za_v;
3168  Numeric lstep;
3169  Index endface;
3170 
3171  do_gridcell_2d_byltest( r_v, lat_v, za_v, lstep, endface, r_start,
3172  lat_start, za_start, -1, 0, ppc, lmax,
3173  lat1, lat3, r1a, r3a, r3b, r1b,
3174  rsurface1, rsurface3 );
3175  // Fill *ppath*
3176  const Index np = r_v.nelem();
3177  ppath_end_2d( ppath, r_v, lat_v, za_v, Vector(np-1,lstep), Vector(np,1),
3178  Vector(np,1), lat_grid, z_field, refellipsoid, ip, ilat,
3179  endface, ppc );
3180 }
3181 
3182 
3183 
3185 
3192  Vector& r_v,
3193  Vector& lat_v,
3194  Vector& lon_v,
3195  Vector& za_v,
3196  Vector& aa_v,
3197  Numeric& lstep,
3198  Index& endface,
3199  const Numeric& r_start0,
3200  const Numeric& lat_start0,
3201  const Numeric& lon_start0,
3202  const Numeric& za_start,
3203  const Numeric& aa_start,
3204  const Numeric& l_start,
3205  const Index& icall,
3206  const Numeric& ppc,
3207  const Numeric& lmax,
3208  const Numeric& lat1,
3209  const Numeric& lat3,
3210  const Numeric& lon5,
3211  const Numeric& lon6,
3212  const Numeric& r15a,
3213  const Numeric& r35a,
3214  const Numeric& r36a,
3215  const Numeric& r16a,
3216  const Numeric& r15b,
3217  const Numeric& r35b,
3218  const Numeric& r36b,
3219  const Numeric& r16b,
3220  const Numeric& rsurface15,
3221  const Numeric& rsurface35,
3222  const Numeric& rsurface36,
3223  const Numeric& rsurface16 )
3224 {
3225  Numeric r_start = r_start0;
3226  Numeric lat_start = lat_start0;
3227  Numeric lon_start = lon_start0;
3228 
3229  assert( icall < 4 );
3230 
3231  // Assert latitude and longitude
3232  assert( lat_start >= lat1 - LATLONTOL );
3233  assert( lat_start <= lat3 + LATLONTOL );
3234  assert( !( abs( lat_start) < POLELAT && lon_start < lon5 - LATLONTOL ) );
3235  assert( !( abs( lat_start) < POLELAT && lon_start > lon6 + LATLONTOL ) );
3236 
3237  // Shift latitude and longitude if outside
3238  if( lat_start < lat1 )
3239  { lat_start = lat1; }
3240  else if( lat_start > lat3 )
3241  { lat_start = lat3; }
3242  if( lon_start < lon5 )
3243  { lon_start = lon5; }
3244  else if( lon_start > lon6 )
3245  { lon_start = lon6; }
3246 
3247  // Radius of lower and upper pressure level at the start position
3248  Numeric rlow = rsurf_at_latlon( lat1, lat3, lon5, lon6,
3249  r15a, r35a, r36a, r16a, lat_start, lon_start );
3250  Numeric rupp = rsurf_at_latlon( lat1, lat3, lon5, lon6,
3251  r15b, r35b, r36b, r16b, lat_start, lon_start );
3252 
3253  // Assert radius (some extra tolerance is needed for radius)
3254  assert( r_start >= rlow - RTOL );
3255  assert( r_start <= rupp + RTOL );
3256 
3257  // Shift radius if outside
3258  if( r_start < rlow )
3259  { r_start = rlow; }
3260  else if( r_start > rupp )
3261  { r_start = rupp; }
3262 
3263  // Position and LOS in cartesian coordinates
3264  Numeric x, y, z, dx, dy, dz;
3265  poslos2cart( x, y, z, dx, dy, dz, r_start, lat_start, lon_start,
3266  za_start, aa_start );
3267 
3268  // Some booleans to check for recursive call
3269  bool unsafe = false;
3270  bool do_surface = false;
3271 
3272  // Determine the position of the end point
3273  //
3274  endface = 0;
3275  //
3276  Numeric r_end, lat_end, lon_end, l_end;
3277 
3278  // Zenith and nadir looking are handled as special cases
3279 
3280  // Zenith looking
3281  if( za_start < ANGTOL )
3282  {
3283  l_end = rupp - r_start;
3284  endface = 4;
3285  }
3286 
3287  // Nadir looking
3288  else if( za_start > 180-ANGTOL )
3289  {
3290  const Numeric rsurface = rsurf_at_latlon( lat1, lat3, lon5, lon6,
3291  rsurface15, rsurface35, rsurface36, rsurface16,
3292  lat_start, lon_start );
3293 
3294  if( rlow > rsurface )
3295  {
3296  l_end = r_start - rlow;
3297  endface = 2;
3298  }
3299  else
3300  {
3301  l_end = r_start - rsurface;
3302  endface = 7;
3303  }
3304  }
3305 
3306  else
3307  {
3308  unsafe = true;
3309 
3310  // Calculate correction terms for the position to compensate for
3311  // numerical inaccuracy.
3312  //
3313  Numeric r_corr, lat_corr, lon_corr;
3314  //
3315  cart2sph( r_corr, lat_corr, lon_corr, x, y, z, lat_start, lon_start,
3316  za_start, aa_start );
3317  //
3318  r_corr -= r_start;
3319  lat_corr -= lat_start;
3320  lon_corr -= lon_start;
3321 
3322  // The end point is found by testing different step lengths until the
3323  // step length has been determined by a precision of *l_acc*.
3324  //
3325  // The first step is to found a length that goes outside the grid cell,
3326  // to find an upper limit. The lower limit is set to 0. The upper
3327  // limit is either set to l_start or it is estimated from the size of
3328  // the grid cell.
3329  // The search algorith is bisection, the next length to test is the
3330  // mean value of the minimum and maximum length limits.
3331  //
3332  if( l_start > 0 )
3333  { l_end = l_start; }
3334  else
3335  { l_end = 2 * ( rupp - rlow ); }
3336  //
3337  Numeric l_in = 0, l_out = l_end;
3338  bool ready = false, startup = true;
3339 
3340  if( rsurface15+RTOL >= r15a || rsurface35+RTOL >= r35a ||
3341  rsurface36+RTOL >= r36a || rsurface16+RTOL >= r16a )
3342  { do_surface = true; }
3343 
3344  while( !ready )
3345  {
3346  cart2sph( r_end, lat_end, lon_end, x+dx*l_end, y+dy*l_end,
3347  z+dz*l_end, lat_start, lon_start, za_start, aa_start );
3348  r_end -= r_corr;
3349  lat_end -= lat_corr;
3350  lon_end -= lon_corr;
3351  resolve_lon( lon_end, lon5, lon6 );
3352 
3353  // Special fix for north-south observations
3354  if( abs( lat_start ) < POLELAT && abs( lat_end ) < POLELAT &&
3355  ( abs(aa_start) < ANGTOL || abs(aa_start) > 180-ANGTOL ) )
3356  { lon_end = lon_start; }
3357 
3358  bool inside = true;
3359 
3360  rlow = rsurf_at_latlon( lat1, lat3, lon5, lon6,
3361  r15a, r35a, r36a, r16a, lat_end, lon_end );
3362 
3363  if( do_surface )
3364  {
3365  const Numeric r_surface =
3366  rsurf_at_latlon( lat1, lat3, lon5, lon6,
3367  rsurface15, rsurface35, rsurface36, rsurface16,
3368  lat_end, lon_end );
3369  if( r_surface >= rlow && r_end <= r_surface )
3370  { inside = false; endface = 7; }
3371  }
3372 
3373  if( inside )
3374  {
3375  if( lat_end < lat1 )
3376  { inside = false; endface = 1; }
3377  else if( lat_end > lat3 )
3378  { inside = false; endface = 3; }
3379  else if( lon_end < lon5 )
3380  { inside = false; endface = 5; }
3381  else if( lon_end > lon6 )
3382  { inside = false; endface = 6; }
3383  else if( r_end < rlow )
3384  { inside = false; endface = 2; }
3385  else
3386  {
3387  rupp = rsurf_at_latlon( lat1, lat3, lon5, lon6,
3388  r15b, r35b, r36b, r16b, lat_end, lon_end );
3389  if( r_end > rupp )
3390  { inside = false; endface = 4; }
3391  }
3392  }
3393 
3394  if( startup )
3395  {
3396  if( inside )
3397  {
3398  l_in = l_end;
3399  l_end *= 5;
3400  }
3401  else
3402  {
3403  l_out = l_end;
3404  l_end = ( l_out + l_in ) / 2;
3405  startup = false;
3406  }
3407  }
3408  else
3409  {
3410  if( inside )
3411  { l_in = l_end; }
3412  else
3413  { l_out = l_end; }
3414 
3415  if( ( l_out - l_in ) < LACC )
3416  { ready = true; }
3417  else
3418  { l_end = ( l_out + l_in ) / 2; }
3419  }
3420  }
3421  }
3422 
3423  //--- Create return vectors
3424  //
3425  Index n = 1;
3426  //
3427  if( lmax > 0 )
3428  {
3429  n = Index( ceil( abs( l_end / lmax ) ) );
3430  if( n < 1 )
3431  { n = 1; }
3432  }
3433  //
3434  r_v.resize( n+1 );
3435  lat_v.resize( n+1 );
3436  lon_v.resize( n+1 );
3437  za_v.resize( n+1 );
3438  aa_v.resize( n+1 );
3439  //
3440  r_v[0] = r_start;
3441  lat_v[0] = lat_start;
3442  lon_v[0] = lon_start;
3443  za_v[0] = za_start;
3444  aa_v[0] = aa_start;
3445  //
3446  lstep = l_end / (Numeric)n;
3447  Numeric l;
3448  bool ready = true;
3449  //
3450  for( Index j=1; j<=n; j++ )
3451  {
3452  l = lstep * (Numeric)j;
3453  cart2poslos( r_v[j], lat_v[j], lon_v[j], za_v[j], aa_v[j], x+dx*l,
3454  y+dy*l, z+dz*l, dx, dy, dz, ppc, lat_start, lon_start,
3455  za_start, aa_start );
3456 
3457  // Shall lon values be shifted (value 0 is already OK)?
3458  resolve_lon( lon_v[j], lon5, lon6 );
3459 
3460  if( j < n )
3461  {
3462  if( unsafe )
3463  {
3464  // Check that r_v[j] is above lower pressure level and the
3465  // surface. This can fail around tangent points. For p-levels
3466  // with constant r this is easy to handle analytically, but the
3467  // problem is tricky in the general case with a non-spherical
3468  // geometry, and this crude solution is used instead. Not the
3469  // most elegant solution, but it works! Added later the same
3470  // check for upper level, after getting assert in that direction.
3471  // The z_field was crazy, but still formerly correct.
3472  rlow = rsurf_at_latlon( lat1, lat3, lon5, lon6, r15a, r35a,
3473  r36a, r16a, lat_v[j], lon_v[j] );
3474  if( do_surface )
3475  {
3476  const Numeric r_surface = rsurf_at_latlon(
3477  lat1, lat3, lon5, lon6, rsurface15, rsurface35,
3478  rsurface36, rsurface16, lat_v[j], lon_v[j] );
3479  const Numeric r_test = max( r_surface, rlow );
3480  if( r_v[j] < r_test )
3481  { ready = false; break; }
3482  }
3483  else if( r_v[j] < rlow )
3484  { ready = false; break; }
3485 
3486  rupp = rsurf_at_latlon( lat1, lat3, lon5, lon6, r15b, r35b,
3487  r36b, r16b, lat_v[j], lon_v[j] );
3488  if( r_v[j] > rupp )
3489  { ready = false; break; }
3490  }
3491  }
3492  else // j==n
3493  {
3494  if( unsafe )
3495  {
3496  // Set end point to be consistent with found endface.
3497  //
3498  if( endface == 1 )
3499  { lat_v[n] = lat1; }
3500  else if( endface == 2 )
3501  { r_v[n] = rsurf_at_latlon( lat1, lat3, lon5, lon6,
3502  r15a, r35a, r36a, r16a,
3503  lat_v[n], lon_v[n] ); }
3504  else if( endface == 3 )
3505  { lat_v[n] = lat3; }
3506  else if( endface == 4 )
3507  { r_v[n] = rsurf_at_latlon( lat1, lat3, lon5, lon6,
3508  r15b, r35b, r36b, r16b,
3509  lat_v[n], lon_v[n] ); }
3510  else if( endface == 5 )
3511  { lon_v[n] = lon5; }
3512  else if( endface == 6 )
3513  { lon_v[n] = lon6; }
3514  else if( endface == 7 )
3515  { r_v[n] = rsurf_at_latlon( lat1, lat3, lon5, lon6,
3516  rsurface15, rsurface35,
3517  rsurface36, rsurface16,
3518  lat_v[n], lon_v[n] ); }
3519  }
3520  }
3521  }
3522 
3523  if( !ready )
3524  { // If an "outside" point found, restart with l as start search length
3525  do_gridcell_3d_byltest( r_v, lat_v, lon_v, za_v, aa_v, lstep, endface,
3526  r_start, lat_start, lon_start, za_start, aa_start,
3527  l, icall+1, ppc, lmax, lat1, lat3, lon5, lon6,
3528  r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b,
3529  rsurface15, rsurface35, rsurface36, rsurface16 );
3530  }
3531 }
3532 
3533 
3534 
3536 
3554  Ppath& ppath,
3555  ConstVectorView lat_grid,
3556  ConstVectorView lon_grid,
3557  ConstTensor3View z_field,
3558  ConstVectorView refellipsoid,
3559  ConstMatrixView z_surface,
3560  const Numeric& lmax )
3561 {
3562  // Radius, zenith angle and latitude of start point.
3563  Numeric r_start, lat_start, lon_start, za_start, aa_start;
3564 
3565  // Lower grid index for the grid cell of interest.
3566  Index ip, ilat, ilon;
3567 
3568  // Radius for corner points, latitude and longitude of the grid cell
3569  //
3570  Numeric lat1, lat3, lon5, lon6;
3571  Numeric r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b;
3572  Numeric rsurface15, rsurface35, rsurface36, rsurface16;
3573 
3574  // Determine the variables defined above and make all possible asserts
3575  ppath_start_3d( r_start, lat_start, lon_start, za_start, aa_start,
3576  ip, ilat, ilon, lat1, lat3, lon5, lon6,
3577  r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b,
3578  rsurface15, rsurface35, rsurface36, rsurface16,
3579  ppath, lat_grid, lon_grid, z_field, refellipsoid, z_surface );
3580 
3581 
3582  // If the field "constant" is negative, this is the first call of the
3583  // function and the path constant shall be calculated.
3584  Numeric ppc;
3585  if( ppath.constant < 0 )
3586  { ppc = geometrical_ppc( r_start, za_start ); }
3587  else
3588  { ppc = ppath.constant; }
3589 
3590 
3591  // Vars to hold found path points, path step length and coding for end face
3592  Vector r_v, lat_v, lon_v, za_v, aa_v;
3593  Numeric lstep;
3594  Index endface;
3595 
3596  do_gridcell_3d_byltest( r_v, lat_v, lon_v, za_v, aa_v, lstep, endface,
3597  r_start, lat_start, lon_start, za_start, aa_start,
3598  -1, 0, ppc, lmax, lat1, lat3, lon5, lon6,
3599  r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b,
3600  rsurface15, rsurface35, rsurface36, rsurface16 );
3601 
3602  // Fill *ppath*
3603  const Index np = r_v.nelem();
3604  ppath_end_3d( ppath, r_v, lat_v, lon_v, za_v, aa_v, Vector(np-1,lstep),
3605  Vector(np,1), Vector(np,1), lat_grid, lon_grid, z_field,
3606  refellipsoid, ip, ilat, ilon, endface, ppc );
3607 }
3608 
3609 
3610 
3611 
3612 
3613 /*===========================================================================
3614  === Core functions for refraction *ppath_step* functions
3615  ===========================================================================*/
3616 
3618 
3659  Workspace& ws,
3660  Array<Numeric>& r_array,
3661  Array<Numeric>& lat_array,
3662  Array<Numeric>& za_array,
3663  Array<Numeric>& l_array,
3664  Array<Numeric>& n_array,
3665  Array<Numeric>& ng_array,
3666  Index& endface,
3667  ConstVectorView p_grid,
3668  ConstVectorView refellipsoid,
3669  ConstTensor3View z_field,
3670  ConstTensor3View t_field,
3671  ConstTensor4View vmr_field,
3672  ConstVectorView f_grid,
3673  const Numeric& lmax,
3674  const Agenda& refr_index_air_agenda,
3675  const Numeric& lraytrace,
3676  const Numeric& rsurface,
3677  const Numeric& r1,
3678  const Numeric& r3,
3679  Numeric r,
3680  Numeric lat,
3681  Numeric za )
3682 {
3683  // Loop boolean
3684  bool ready = false;
3685 
3686  // Store first point
3687  Numeric refr_index_air, refr_index_air_group;
3688  get_refr_index_1d( ws, refr_index_air, refr_index_air_group,
3689  refr_index_air_agenda, p_grid, refellipsoid,
3690  z_field, t_field, vmr_field, f_grid, r );
3691  r_array.push_back( r );
3692  lat_array.push_back( lat );
3693  za_array.push_back( za );
3694  n_array.push_back( refr_index_air );
3695  ng_array.push_back( refr_index_air_group );
3696 
3697  // Variables for output from do_gridcell_2d
3698  Vector r_v, lat_v, za_v;
3699  Numeric lstep, lcum = 0, dlat;
3700 
3701  while( !ready )
3702  {
3703  // Constant for the geometrical step to make
3704  const Numeric ppc_step = geometrical_ppc( r, za );
3705 
3706  // Where will a geometric path exit the grid cell?
3707  do_gridrange_1d( r_v, lat_v, za_v, lstep, endface, r, lat, za,
3708  ppc_step, -1, r1, r3, rsurface );
3709  assert( r_v.nelem() == 2 );
3710 
3711  // If *lstep* is <= *lraytrace*, extract the found end point (if not
3712  // a tangent point, we are ready).
3713  // Otherwise, we make a geometrical step with length *lraytrace*.
3714 
3715  Numeric za_flagside = za;
3716 
3717  if( lstep <= lraytrace )
3718  {
3719  r = r_v[1];
3720  dlat = lat_v[1] - lat;
3721  lat = lat_v[1];
3722  lcum += lstep;
3723  ready = true;
3724  }
3725  else
3726  {
3727  Numeric l;
3728  if( za <= 90 )
3729  { l = geompath_l_at_r(ppc_step,r) + lraytrace; }
3730  else
3731  {
3732  l = geompath_l_at_r(ppc_step,r) - lraytrace;
3733  if( l < 0 )
3734  { za_flagside = 180-za_flagside; } // Tangent point passed!
3735  }
3736 
3737  r = geompath_r_at_l( ppc_step, l );
3738 
3739  const Numeric lat_new = geompath_lat_at_za( za, lat,
3740  geompath_za_at_r( ppc_step, za_flagside, r) );
3741  dlat = lat_new - lat;
3742  lat = lat_new;
3743  lstep = lraytrace;
3744  lcum += lraytrace;
3745  }
3746 
3747  // Refractive index at new point
3748  Numeric dndr;
3749  refr_gradients_1d( ws, refr_index_air, refr_index_air_group, dndr,
3750  refr_index_air_agenda, p_grid, refellipsoid, z_field,
3751  t_field, vmr_field, f_grid, r );
3752 
3753  // Calculate LOS zenith angle at found point.
3754  const Numeric za_rad = DEG2RAD * za;
3755  //
3756  za += - dlat; // Pure geometrical effect
3757  //
3758  za += (RAD2DEG * lstep / refr_index_air) * ( -sin(za_rad) * dndr );
3759 
3760  // Make sure that obtained *za* is inside valid range
3761  if( za < 0 )
3762  { za = -za; }
3763  else if( za > 180 )
3764  { za -= 360; }
3765 
3766  // Store found point?
3767  if( ready || ( lmax > 0 && lcum + lraytrace > lmax ) )
3768  {
3769  r_array.push_back( r );
3770  lat_array.push_back( lat );
3771  za_array.push_back( za );
3772  n_array.push_back( refr_index_air );
3773  ng_array.push_back( refr_index_air_group );
3774  l_array.push_back( lcum );
3775  lcum = 0;
3776  }
3777  }
3778 }
3779 
3780 
3782 
3810  Workspace& ws,
3811  Ppath& ppath,
3812  ConstVectorView p_grid,
3813  ConstTensor3View z_field,
3814  ConstTensor3View t_field,
3815  ConstTensor4View vmr_field,
3816  ConstVectorView f_grid,
3817  ConstVectorView refellipsoid,
3818  const Numeric& z_surface,
3819  const Numeric& lmax,
3820  const Agenda& refr_index_air_agenda,
3821  const String& rtrace_method,
3822  const Numeric& lraytrace )
3823 {
3824  // Starting radius, zenith angle and latitude
3825  Numeric r_start, lat_start, za_start;
3826 
3827  // Index of the pressure level being the lower limit for the
3828  // grid range of interest.
3829  Index ip;
3830 
3831  // Determine the variables defined above, and make asserts of input
3832  ppath_start_1d( r_start, lat_start, za_start, ip, ppath );
3833 
3834  // If the field "constant" is negative, this is the first call of the
3835  // function and the path constant shall be calculated.
3836  // If the sensor is placed outside the atmosphere, the constant is
3837  // already set.
3838  Numeric ppc;
3839  if( ppath.constant < 0 )
3840  {
3841  Numeric refr_index_air, refr_index_air_group;
3842  get_refr_index_1d( ws, refr_index_air, refr_index_air_group,
3843  refr_index_air_agenda, p_grid, refellipsoid,
3844  z_field, t_field, vmr_field, f_grid, r_start );
3845  ppc = refraction_ppc( r_start, za_start, refr_index_air );
3846  }
3847  else
3848  { ppc = ppath.constant; }
3849 
3850 
3851  // Perform the ray tracing
3852  //
3853  // Arrays to store found ray tracing points
3854  // (Vectors don't work here as we don't know how many points there will be)
3855  Array<Numeric> r_array, lat_array, za_array, l_array, n_array, ng_array;
3856  Index endface;
3857  //
3858  if( rtrace_method == "linear_basic" )
3859  {
3860  /*
3861  raytrace_1d_linear_basic( ws, r_array, lat_array, za_array, l_array,
3862  n_array, ng_array, endface, refellipsoid, p_grid, z_field, t_field,
3863  vmr_field, f_grid, lmax, refr_index_air_agenda,
3864  lraytrace, ppc, refellipsoid[0] + z_surface,
3865  refellipsoid[0]+z_field[ip], refellipsoid[0]+z_field[ip+1],
3866  r_start, lat_start, za_start );
3867  */
3868  raytrace_1d_linear_basic( ws, r_array, lat_array, za_array, l_array,
3869  n_array, ng_array, endface, p_grid, refellipsoid, z_field, t_field,
3870  vmr_field, f_grid, lmax, refr_index_air_agenda, lraytrace,
3871  refellipsoid[0] + z_surface, refellipsoid[0]+z_field(ip,0,0),
3872  refellipsoid[0]+z_field(ip+1,0,0), r_start, lat_start, za_start );
3873  }
3874  else
3875  {
3876  // Make sure we fail if called with an invalid rtrace_method.
3877  assert(false);
3878  }
3879 
3880  // Fill *ppath*
3881  //
3882  const Index np = r_array.nelem();
3883  Vector r_v(np), lat_v(np), za_v(np), l_v(np-1), n_v(np), ng_v(np);
3884  for( Index i=0; i<np; i++ )
3885  {
3886  r_v[i] = r_array[i];
3887  lat_v[i] = lat_array[i];
3888  za_v[i] = za_array[i];
3889  n_v[i] = n_array[i];
3890  ng_v[i] = ng_array[i];
3891  if( i < np-1 )
3892  { l_v[i] = l_array[i]; }
3893  }
3894  //
3895  ppath_end_1d( ppath, r_v, lat_v, za_v, l_v, n_v, ng_v, z_field(joker,0,0),
3896  refellipsoid, ip, endface, ppc );
3897 }
3898 
3899 
3900 
3902 
3948  Workspace& ws,
3949  Array<Numeric>& r_array,
3950  Array<Numeric>& lat_array,
3951  Array<Numeric>& za_array,
3952  Array<Numeric>& l_array,
3953  Array<Numeric>& n_array,
3954  Array<Numeric>& ng_array,
3955  Index& endface,
3956  ConstVectorView p_grid,
3957  ConstVectorView lat_grid,
3958  ConstVectorView refellipsoid,
3959  ConstTensor3View z_field,
3960  ConstTensor3View t_field,
3961  ConstTensor4View vmr_field,
3962  ConstVectorView f_grid,
3963  const Numeric& lmax,
3964  const Agenda& refr_index_air_agenda,
3965  const Numeric& lraytrace,
3966  const Numeric& lat1,
3967  const Numeric& lat3,
3968  const Numeric& rsurface1,
3969  const Numeric& rsurface3,
3970  const Numeric& r1a,
3971  const Numeric& r3a,
3972  const Numeric& r3b,
3973  const Numeric& r1b,
3974  Numeric r,
3975  Numeric lat,
3976  Numeric za )
3977 {
3978  // Loop boolean
3979  bool ready = false;
3980 
3981  // Store first point
3982  Numeric refr_index_air, refr_index_air_group;
3983  get_refr_index_2d( ws, refr_index_air, refr_index_air_group,
3984  refr_index_air_agenda, p_grid, lat_grid, refellipsoid,
3985  z_field, t_field, vmr_field, f_grid, r, lat );
3986  r_array.push_back( r );
3987  lat_array.push_back( lat );
3988  za_array.push_back( za );
3989  n_array.push_back( refr_index_air );
3990  ng_array.push_back( refr_index_air_group );
3991 
3992  // Variables for output from do_gridcell_2d
3993  Vector r_v, lat_v, za_v;
3994  Numeric lstep, lcum = 0, dlat;
3995 
3996  while( !ready )
3997  {
3998  // Constant for the geometrical step to make
3999  const Numeric ppc_step = geometrical_ppc( r, za );
4000 
4001  // Where will a geometric path exit the grid cell?
4002  do_gridcell_2d_byltest( r_v, lat_v, za_v, lstep, endface, r, lat, za,
4003  lraytrace, 0, ppc_step, -1, lat1, lat3,
4004  r1a, r3a, r3b, r1b, rsurface1, rsurface3 );
4005  assert( r_v.nelem() == 2 );
4006 
4007  // If *lstep* is <= *lraytrace*, extract the found end point (if not
4008  // a tangent point, we are ready).
4009  // Otherwise, we make a geometrical step with length *lraytrace*.
4010 
4011  Numeric za_flagside = za;
4012 
4013  if( lstep <= lraytrace )
4014  {
4015  r = r_v[1];
4016  dlat = lat_v[1] - lat;
4017  lat = lat_v[1];
4018  lcum += lstep;
4019  ready = true;
4020  }
4021  else
4022  {
4023  Numeric l;
4024  if( abs(za) <= 90 )
4025  { l = geompath_l_at_r(ppc_step,r) + lraytrace; }
4026  else
4027  {
4028  l = geompath_l_at_r(ppc_step,r) - lraytrace;
4029  if( l < 0 ) // Tangent point passed!
4030  { za_flagside = sign(za)*180-za_flagside; }
4031  }
4032 
4033  r = geompath_r_at_l( ppc_step, l );
4034 
4035  const Numeric lat_new = geompath_lat_at_za( za, lat,
4036  geompath_za_at_r( ppc_step, za_flagside, r) );
4037  dlat = lat_new - lat;
4038  lat = lat_new;
4039  lstep = lraytrace;
4040  lcum += lraytrace;
4041 
4042  // For paths along the latitude end faces we can end up outside the
4043  // grid cell. We simply look for points outisde the grid cell.
4044  if( lat < lat1 )
4045  { lat = lat1; }
4046  else if( lat > lat3 )
4047  { lat = lat3; }
4048  }
4049 
4050  // Refractive index at new point
4051  Numeric dndr, dndlat;
4052  refr_gradients_2d( ws, refr_index_air, refr_index_air_group, dndr,
4053  dndlat, refr_index_air_agenda, p_grid, lat_grid,
4054  refellipsoid, z_field, t_field, vmr_field, f_grid,
4055  r, lat );
4056 
4057  // Calculate LOS zenith angle at found point.
4058  const Numeric za_rad = DEG2RAD * za;
4059  //
4060  za += - dlat; // Pure geometrical effect
4061  //
4062  za += (RAD2DEG * lstep / refr_index_air) * ( -sin(za_rad) * dndr +
4063  cos(za_rad) * dndlat );
4064 
4065  // Make sure that obtained *za* is inside valid range
4066  if( za < -180 )
4067  { za += 360; }
4068  else if( za > 180 )
4069  { za -= 360; }
4070 
4071  // If the path is zenith/nadir along a latitude end face, we must check
4072  // that the path does not exit with new *za*.
4073  if( lat == lat1 && za < 0 )
4074  { endface = 1; ready = 1; }
4075  else if( lat == lat3 && za > 0 )
4076  { endface = 3; ready = 1; }
4077 
4078  // Store found point?
4079  if( ready || ( lmax > 0 && lcum + lraytrace > lmax ) )
4080  {
4081  r_array.push_back( r );
4082  lat_array.push_back( lat );
4083  za_array.push_back( za );
4084  n_array.push_back( refr_index_air );
4085  ng_array.push_back( refr_index_air_group );
4086  l_array.push_back( lcum );
4087  lcum = 0;
4088  }
4089  }
4090 }
4091 
4092 
4093 
4095 
4122  Workspace& ws,
4123  Ppath& ppath,
4124  ConstVectorView p_grid,
4125  ConstVectorView lat_grid,
4126  ConstTensor3View z_field,
4127  ConstTensor3View t_field,
4128  ConstTensor4View vmr_field,
4129  ConstVectorView f_grid,
4130  ConstVectorView refellipsoid,
4131  ConstVectorView z_surface,
4132  const Numeric& lmax,
4133  const Agenda& refr_index_air_agenda,
4134  const String& rtrace_method,
4135  const Numeric& lraytrace )
4136 {
4137  // Radius, zenith angle and latitude of start point.
4138  Numeric r_start, lat_start, za_start;
4139 
4140  // Lower grid index for the grid cell of interest.
4141  Index ip, ilat;
4142 
4143  // Radii and latitudes set by *ppath_start_2d*.
4144  Numeric lat1, lat3, r1a, r3a, r3b, r1b, rsurface1, rsurface3;
4145 
4146  // Determine the variables defined above and make all possible asserts
4147  ppath_start_2d( r_start, lat_start, za_start, ip, ilat,
4148  lat1, lat3, r1a, r3a, r3b, r1b, rsurface1, rsurface3,
4149  ppath, lat_grid, z_field(joker,joker,0), refellipsoid,
4150  z_surface );
4151 
4152  // Perform the ray tracing
4153  //
4154  // No constant for the path is valid here.
4155  //
4156  // Arrays to store found ray tracing points
4157  // (Vectors don't work here as we don't know how many points there will be)
4158  Array<Numeric> r_array, lat_array, za_array, l_array, n_array, ng_array;
4159  Index endface;
4160  //
4161  if( rtrace_method == "linear_basic" )
4162  {
4163  raytrace_2d_linear_basic( ws, r_array, lat_array, za_array, l_array,
4164  n_array, ng_array, endface, p_grid, lat_grid,
4165  refellipsoid, z_field, t_field, vmr_field,
4166  f_grid, lmax,
4167  refr_index_air_agenda, lraytrace, lat1, lat3,
4168  rsurface1, rsurface3, r1a, r3a, r3b, r1b,
4169  r_start, lat_start, za_start );
4170  }
4171  else
4172  {
4173  // Make sure we fail if called with an invalid rtrace_method.
4174  assert(false);
4175  }
4176 
4177  // Fill *ppath*
4178  //
4179  const Index np = r_array.nelem();
4180  Vector r_v(np), lat_v(np), za_v(np), l_v(np-1), n_v(np), ng_v(np);
4181  for( Index i=0; i<np; i++ )
4182  {
4183  r_v[i] = r_array[i];
4184  lat_v[i] = lat_array[i];
4185  za_v[i] = za_array[i];
4186  n_v[i] = n_array[i];
4187  ng_v[i] = ng_array[i];
4188  if( i < np-1 )
4189  { l_v[i] = l_array[i]; }
4190  }
4191  //
4192  ppath_end_2d( ppath, r_v, lat_v, za_v, l_v, n_v, ng_v, lat_grid,
4193  z_field(joker,joker,0), refellipsoid, ip, ilat, endface, -1 );
4194 }
4195 
4196 
4197 
4199 
4259  Workspace& ws,
4260  Array<Numeric>& r_array,
4261  Array<Numeric>& lat_array,
4262  Array<Numeric>& lon_array,
4263  Array<Numeric>& za_array,
4264  Array<Numeric>& aa_array,
4265  Array<Numeric>& l_array,
4266  Array<Numeric>& n_array,
4267  Array<Numeric>& ng_array,
4268  Index& endface,
4269  ConstVectorView refellipsoid,
4270  ConstVectorView p_grid,
4271  ConstVectorView lat_grid,
4272  ConstVectorView lon_grid,
4273  ConstTensor3View z_field,
4274  ConstTensor3View t_field,
4275  ConstTensor4View vmr_field,
4276  ConstVectorView f_grid,
4277  const Numeric& lmax,
4278  const Agenda& refr_index_air_agenda,
4279  const Numeric& lraytrace,
4280  const Numeric& lat1,
4281  const Numeric& lat3,
4282  const Numeric& lon5,
4283  const Numeric& lon6,
4284  const Numeric& rsurface15,
4285  const Numeric& rsurface35,
4286  const Numeric& rsurface36,
4287  const Numeric& rsurface16,
4288  const Numeric& r15a,
4289  const Numeric& r35a,
4290  const Numeric& r36a,
4291  const Numeric& r16a,
4292  const Numeric& r15b,
4293  const Numeric& r35b,
4294  const Numeric& r36b,
4295  const Numeric& r16b,
4296  Numeric r,
4297  Numeric lat,
4298  Numeric lon,
4299  Numeric za,
4300  Numeric aa )
4301 {
4302  // Loop boolean
4303  bool ready = false;
4304 
4305  // Store first point
4306  Numeric refr_index_air, refr_index_air_group;
4307  get_refr_index_3d( ws, refr_index_air, refr_index_air_group,
4308  refr_index_air_agenda, p_grid, lat_grid, lon_grid,
4309  refellipsoid, z_field, t_field, vmr_field, f_grid,
4310  r, lat, lon );
4311  r_array.push_back( r );
4312  lat_array.push_back( lat );
4313  lon_array.push_back( lon );
4314  za_array.push_back( za );
4315  aa_array.push_back( aa );
4316  n_array.push_back( refr_index_air );
4317  ng_array.push_back( refr_index_air_group );
4318 
4319  // Variables for output from do_gridcell_3d
4320  Vector r_v, lat_v, lon_v, za_v, aa_v;
4321  Numeric lstep, lcum = 0;
4322  Numeric za_new, aa_new;
4323 
4324  while( !ready )
4325  {
4326  // Constant for the geometrical step to make
4327  const Numeric ppc_step = geometrical_ppc( r, za );
4328 
4329  // Where will a geometric path exit the grid cell?
4330  do_gridcell_3d_byltest( r_v, lat_v, lon_v, za_v, aa_v, lstep, endface,
4331  r, lat, lon, za, aa, lraytrace, 0,
4332  ppc_step, -1, lat1, lat3, lon5, lon6,
4333  r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b,
4334  rsurface15, rsurface35, rsurface36, rsurface16 );
4335  assert( r_v.nelem() == 2 );
4336 
4337  // If *lstep* is <= *lraytrace*, extract the found end point.
4338  // Otherwise, we make a geometrical step with length *lraytrace*.
4339 
4340  if( lstep <= lraytrace )
4341  {
4342  r = r_v[1];
4343  lat = lat_v[1];
4344  lon = lon_v[1];
4345  za_new = za_v[1];
4346  aa_new = aa_v[1];
4347  lcum += lstep;
4348  ready = true;
4349  }
4350  else
4351  {
4352  // Sensor pos and LOS in cartesian coordinates
4353  Numeric x, y, z, dx, dy, dz, lat_new, lon_new;
4354  //
4355  poslos2cart( x, y, z, dx, dy, dz, r, lat, lon, za, aa );
4356  lstep = lraytrace;
4357  cart2poslos( r, lat_new, lon_new, za_new, aa_new, x+dx*lstep,
4358  y+dy*lstep, z+dz*lstep, dx, dy, dz, ppc_step,
4359  lat, lon, za, aa );
4360  lcum += lstep;
4361 
4362  // Shall lon values be shifted?
4363  resolve_lon( lon_new, lon5, lon6 );
4364 
4365  lat = lat_new;
4366  lon = lon_new;
4367  }
4368 
4369  // Refractive index at new point
4370  Numeric dndr, dndlat, dndlon;
4371  refr_gradients_3d( ws, refr_index_air, refr_index_air_group,
4372  dndr, dndlat, dndlon, refr_index_air_agenda,
4373  p_grid, lat_grid, lon_grid, refellipsoid,
4374  z_field, t_field, vmr_field,
4375  f_grid, r, lat, lon );
4376 
4377  // Calculate LOS zenith angle at found point.
4378  const Numeric aterm = RAD2DEG * lstep / refr_index_air;
4379  const Numeric za_rad = DEG2RAD * za;
4380  const Numeric aa_rad = DEG2RAD * aa;
4381  const Numeric sinza = sin( za_rad );
4382  const Numeric sinaa = sin( aa_rad );
4383  const Numeric cosaa = cos( aa_rad );
4384  //
4385  Vector los(2); los[0] = za_new; los[1] = aa_new;
4386  //
4387  if( za < ANGTOL || za > 180-ANGTOL )
4388  {
4389  los[0] += aterm * ( cos(za_rad) *
4390  ( cosaa * dndlat + sinaa * dndlon ) );
4391  los[1] = RAD2DEG * atan2( dndlon, dndlat);
4392  }
4393  else
4394  {
4395  los[0] += aterm * ( -sinza * dndr + cos(za_rad) *
4396  ( cosaa * dndlat + sinaa * dndlon ) );
4397  los[1] += aterm * sinza * ( cosaa * dndlon - sinaa * dndlat );
4398  }
4399  //
4400  adjust_los( los, 3 );
4401  //
4402  za = los[0];
4403  aa = los[1];
4404 
4405  // For some cases where the path goes along an end face,
4406  // it could be the case that the refraction bends the path out
4407  // of the grid cell.
4408  if( za > 0 && za < 180 )
4409  {
4410  if( lon == lon5 && aa < 0 )
4411  { endface = 5; ready = 1; }
4412  else if( lon == lon6 && aa > 0 )
4413  { endface = 6; ready = 1; }
4414  else if( lat == lat1 && lat != -90 && abs( aa ) > 90 )
4415  { endface = 1; ready = 1; }
4416  else if( lat == lat3 && lat != 90 && abs( aa ) < 90 )
4417  { endface = 3; ready = 1; }
4418  }
4419 
4420  // Store found point?
4421  if( ready || ( lmax > 0 && lcum + lraytrace > lmax ) )
4422  {
4423  r_array.push_back( r );
4424  lat_array.push_back( lat );
4425  lon_array.push_back( lon );
4426  za_array.push_back( za );
4427  aa_array.push_back( aa );
4428  n_array.push_back( refr_index_air );
4429  ng_array.push_back( refr_index_air_group );
4430  l_array.push_back( lcum );
4431  lcum = 0;
4432  }
4433  }
4434 }
4435 
4436 
4437 
4439 
4467  Workspace& ws,
4468  Ppath& ppath,
4469  ConstVectorView p_grid,
4470  ConstVectorView lat_grid,
4471  ConstVectorView lon_grid,
4472  ConstTensor3View z_field,
4473  ConstTensor3View t_field,
4474  ConstTensor4View vmr_field,
4475  ConstVectorView f_grid,
4476  ConstVectorView refellipsoid,
4477  ConstMatrixView z_surface,
4478  const Numeric& lmax,
4479  const Agenda& refr_index_air_agenda,
4480  const String& rtrace_method,
4481  const Numeric& lraytrace )
4482 {
4483  // Radius, zenith angle and latitude of start point.
4484  Numeric r_start, lat_start, lon_start, za_start, aa_start;
4485 
4486  // Lower grid index for the grid cell of interest.
4487  Index ip, ilat, ilon;
4488 
4489  // Radius for corner points, latitude and longitude of the grid cell
4490  //
4491  Numeric lat1, lat3, lon5, lon6;
4492  Numeric r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b;
4493  Numeric rsurface15, rsurface35, rsurface36, rsurface16;
4494 
4495  // Determine the variables defined above and make all possible asserts
4496  ppath_start_3d( r_start, lat_start, lon_start, za_start, aa_start,
4497  ip, ilat, ilon, lat1, lat3, lon5, lon6,
4498  r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b,
4499  rsurface15, rsurface35, rsurface36, rsurface16,
4500  ppath, lat_grid, lon_grid, z_field, refellipsoid, z_surface );
4501 
4502  // Perform the ray tracing
4503  //
4504  // No constant for the path is valid here.
4505  //
4506  // Arrays to store found ray tracing points
4507  // (Vectors don't work here as we don't know how many points there will be)
4508  Array<Numeric> r_array, lat_array, lon_array, za_array, aa_array;
4509  Array<Numeric> l_array, n_array, ng_array;
4510  Index endface;
4511  //
4512  if( rtrace_method == "linear_basic" )
4513  {
4514  raytrace_3d_linear_basic( ws, r_array, lat_array, lon_array, za_array,
4515  aa_array, l_array, n_array, ng_array, endface,
4516  refellipsoid, p_grid, lat_grid, lon_grid,
4517  z_field, t_field, vmr_field,
4518  f_grid, lmax, refr_index_air_agenda, lraytrace,
4519  lat1, lat3, lon5, lon6,
4520  rsurface15, rsurface35, rsurface36, rsurface16,
4521  r15a, r35a, r36a, r16a, r15b, r35b, r36b, r16b,
4522  r_start, lat_start, lon_start, za_start, aa_start );
4523  }
4524  else
4525  {
4526  // Make sure we fail if called with an invalid rtrace_method.
4527  assert(false);
4528  }
4529 
4530  // Fill *ppath*
4531  //
4532  const Index np = r_array.nelem();
4533  Vector r_v(np), lat_v(np), lon_v(np), za_v(np), aa_v(np), l_v(np-1);
4534  Vector n_v(np), ng_v(np);
4535  for( Index i=0; i<np; i++ )
4536  {
4537  r_v[i] = r_array[i];
4538  lat_v[i] = lat_array[i];
4539  lon_v[i] = lon_array[i];
4540  za_v[i] = za_array[i];
4541  aa_v[i] = aa_array[i];
4542  n_v[i] = n_array[i];
4543  ng_v[i] = ng_array[i];
4544  if( i < np-1 )
4545  { l_v[i] = l_array[i]; }
4546  }
4547  //
4548  // Fill *ppath*
4549  ppath_end_3d( ppath, r_v, lat_v, lon_v, za_v, aa_v, l_v, n_v, ng_v, lat_grid,
4550  lon_grid, z_field, refellipsoid, ip, ilat, ilon, endface, -1 );
4551 }
4552 
4553 
4554 
4555 
4556 
4557 /*===========================================================================
4558  === Main functions
4559  ===========================================================================*/
4560 
4562 
4607  Ppath& ppath,
4608  const Index& atmosphere_dim,
4609  ConstVectorView p_grid,
4610  ConstVectorView lat_grid,
4611  ConstVectorView lon_grid,
4612  ConstTensor3View z_field,
4613  ConstVectorView refellipsoid,
4614  ConstMatrixView z_surface,
4615  const Index& cloudbox_on,
4616  const ArrayOfIndex& cloudbox_limits,
4617  const bool& ppath_inside_cloudbox_do,
4618  ConstVectorView rte_pos,
4619  ConstVectorView rte_los,
4620  const Verbosity& verbosity)
4621 {
4622  CREATE_OUT1;
4623 
4624  // This function contains no checks or asserts as it is only a sub-function.
4625 
4626  // Allocate the ppath structure
4627  ppath_init_structure( ppath, atmosphere_dim, 1 );
4628 
4629  // Index of last pressure level
4630  const Index lp = p_grid.nelem() - 1;
4631 
4632  // The different atmospheric dimensionalities are handled seperately
4633 
4634  //-- 1D ---------------------------------------------------------------------
4635  if( atmosphere_dim == 1 )
4636  {
4637  // End position and LOS
4638  ppath.end_pos[0] = rte_pos[0];
4639  ppath.end_pos[1] = 0;
4640  ppath.end_los[0] = rte_los[0];
4641 
4642  // Sensor is inside the model atmosphere:
4643  if( rte_pos[0] < z_field(lp,0,0) )
4644  {
4645  // Check that the sensor is above the surface
4646  if( (rte_pos[0] + RTOL) < z_surface(0,0) )
4647  {
4648  ostringstream os;
4649  os << "The ppath starting point is placed "
4650  << (z_surface(0,0) - rte_pos[0])/1e3
4651  << " km below the surface.";
4652  throw runtime_error(os.str());
4653  }
4654 
4655  // Set ppath
4656  ppath.pos(0,joker) = ppath.end_pos;
4657  ppath.r[0] = refellipsoid[0] + rte_pos[0];
4658  ppath.los(0,joker) = ppath.end_los;
4659  //
4660  gridpos( ppath.gp_p, z_field(joker,0,0), ppath.pos(0,0) );
4661  gridpos_check_fd( ppath.gp_p[0] );
4662 
4663  // Is the sensor on the surface looking down?
4664  // If yes and the sensor is inside the cloudbox, the background will
4665  // be changed below.
4666  if( ppath.pos(0,0) <= z_surface(0,0) && ppath.los(0,0) > 90 )
4667  { ppath_set_background( ppath, 2 ); }
4668 
4669  // Outside cloudbox:
4670  // Check sensor position with respect to cloud box.
4671  if( cloudbox_on && !ppath_inside_cloudbox_do )
4672  {
4673  const Numeric fgp = fractional_gp( ppath.gp_p[0] );
4674  if( fgp >= cloudbox_limits[0] && fgp <= cloudbox_limits[1] )
4675  { ppath_set_background( ppath, 4 ); }
4676  }
4677 
4678  // Inside cloudbox:
4679  DEBUG_ONLY( if( ppath_inside_cloudbox_do )
4680  {
4681  const Numeric fgp = fractional_gp( ppath.gp_p[0] );
4682  assert( fgp >= cloudbox_limits[0] && fgp <= cloudbox_limits[1] );
4683  } )
4684  }
4685 
4686  // Sensor is outside the model atmosphere:
4687  else
4688  {
4689  // We can here set ppc and n as we are outside the atmosphere
4690  ppath.nreal = 1.0;
4691  ppath.ngroup = 1.0;
4692  ppath.constant = geometrical_ppc( refellipsoid[0] + rte_pos[0],
4693  rte_los[0] );
4694 
4695  // Totally outside
4696  if( rte_los[0] <= 90 ||
4697  ppath.constant >= refellipsoid[0] + z_field(lp,0,0) )
4698  {
4699  ppath.pos(0,0) = rte_pos[0];
4700  ppath.pos(0,1) = 0;
4701  ppath.r[0] = refellipsoid[0] + rte_pos[0];
4702  ppath.los(0,0) = rte_los[0];
4703  //
4704  ppath_set_background( ppath, 1 );
4705  out1 << " --- WARNING ---, path is totally outside of the "
4706  << "model atmosphere\n";
4707  }
4708 
4709  // Path enters the atmosphere.
4710  else
4711  {
4712  ppath.r[0] = refellipsoid[0] + z_field(lp,0,0);
4713  ppath.pos(0,0) = z_field(lp,0,0);
4714  ppath.los(0,0) = geompath_za_at_r( ppath.constant, rte_los[0],
4715  ppath.r[0] );
4716  ppath.pos(0,1) = geompath_lat_at_za( rte_los[0], 0,
4717  ppath.los(0,0) );
4718  ppath.end_lstep = geompath_l_at_r( ppath.constant,
4719  refellipsoid[0] + rte_pos[0] ) -
4720  geompath_l_at_r( ppath.constant, ppath.r[0] );
4721 
4722  // Here we know the grid position exactly
4723  ppath.gp_p[0].idx = lp-1;
4724  ppath.gp_p[0].fd[0] = 1;
4725  ppath.gp_p[0].fd[1] = 0;
4726 
4727  // If cloud box reaching TOA, we have also found the background
4728  if( cloudbox_on && cloudbox_limits[1] == lp )
4729  { ppath_set_background( ppath, 3 ); }
4730  }
4731  }
4732  } // End 1D
4733 
4734 
4735  //-- 2D ---------------------------------------------------------------------
4736  else if( atmosphere_dim == 2 )
4737  {
4738  // End position and LOS
4739  ppath.end_pos[0] = rte_pos[0];
4740  ppath.end_pos[1] = rte_pos[1];
4741  ppath.end_los[0] = rte_los[0];
4742 
4743  // Index of last latitude
4744  const Index llat = lat_grid.nelem() -1;
4745 
4746  // Is sensor inside range of lat_grid?
4747  // If yes, determine TOA altitude at sensor position
4748  GridPos gp_lat;
4749  Vector itw(2);
4750  bool islatin = false;
4751  Numeric r_e; // Ellipsoid radius at sensor position
4752  Numeric z_toa = -99e99;
4753  if( rte_pos[1] > lat_grid[0] && rte_pos[1] < lat_grid[llat] )
4754  {
4755  islatin = true;
4756  gridpos( gp_lat, lat_grid, rte_pos[1] );
4757  interpweights( itw, gp_lat );
4758  z_toa = interp( itw, z_field(lp,joker,0), gp_lat );
4759  r_e = refell2d( refellipsoid, lat_grid, gp_lat );
4760  }
4761  else
4762  { r_e = refell2r( refellipsoid, rte_pos[1] ); }
4763 
4764  // Sensor is inside the model atmosphere:
4765  if( islatin && rte_pos[0] < z_toa )
4766  {
4767  const Numeric z_s = interp( itw, z_surface(joker,0), gp_lat );
4768 
4769  // Check that the sensor is above the surface
4770  if( (rte_pos[0] + RTOL) < z_s )
4771  {
4772  ostringstream os;
4773  os << "The ppath starting point is placed "
4774  << (z_s - rte_pos[0])/1e3 << " km below the surface.";
4775  throw runtime_error(os.str());
4776  }
4777 
4778  // Set ppath
4779  ppath.pos(0,joker) = ppath.end_pos;
4780  ppath.r[0] = r_e + rte_pos[0];
4781  ppath.los(0,joker) = ppath.end_los;
4782 
4783  // Determine gp_p (gp_lat is recalculated, but ...)
4784  GridPos gp_lon_dummy;
4785  rte_pos2gridpos( ppath.gp_p[0], ppath.gp_lat[0], gp_lon_dummy,
4786  atmosphere_dim, p_grid, lat_grid, lon_grid, z_field, rte_pos );
4787  gridpos_check_fd( ppath.gp_p[0] );
4788  gridpos_check_fd( ppath.gp_lat[0] );
4789 
4790  // Is the sensor on the surface looking down?
4791  // If yes and the sensor is inside the cloudbox, the background will
4792  // be changed below.
4793  if( ppath.pos(0,0) <= z_s )
4794  {
4795  // Calculate radial slope of the surface
4796  Numeric rslope;
4797  plevel_slope_2d( rslope, lat_grid, refellipsoid,
4798  z_surface(joker,0), gp_lat, ppath.los(0,0) );
4799 
4800  // Calculate angular tilt of the surface
4801  const Numeric atilt = plevel_angletilt( r_e + z_s, rslope );
4802 
4803  // Are we looking down into the surface?
4804  // If yes and the sensor is inside the cloudbox, the background
4805  // will be changed below.
4806  if( is_los_downwards( ppath.los(0,0), atilt ) )
4807  { ppath_set_background( ppath, 2 ); }
4808  }
4809 
4810  // Outside cloudbox:
4811  // Check sensor position with respect to cloud box.
4812  if( cloudbox_on && !ppath_inside_cloudbox_do )
4813  {
4814  const Numeric fgp = fractional_gp( ppath.gp_p[0] );
4815  const Numeric fgl = fractional_gp( ppath.gp_lat[0] );
4816  if( fgp >= cloudbox_limits[0] && fgp <= cloudbox_limits[1] &&
4817  fgl >= cloudbox_limits[2] && fgl <= cloudbox_limits[3] )
4818  { ppath_set_background( ppath, 4 ); }
4819  }
4820 
4821  // Inside cloudbox:
4822  DEBUG_ONLY( if( ppath_inside_cloudbox_do )
4823  {
4824  const Numeric fgp = fractional_gp( ppath.gp_p[0] );
4825  const Numeric fgl = fractional_gp( ppath.gp_lat[0] );
4826  assert( fgp >= cloudbox_limits[0] && fgp <= cloudbox_limits[1] &&
4827  fgl >= cloudbox_limits[2] && fgl <= cloudbox_limits[3] );
4828  } )
4829  }
4830 
4831  // Sensor is outside the model atmosphere:
4832  else
4833  {
4834  // Handle cases when the sensor looks in the wrong way
4835  if( ( rte_pos[1] <= lat_grid[0] && rte_los[0] <= 0 ) ||
4836  ( rte_pos[1] >= lat_grid[llat] && rte_los[0] >= 0 ) )
4837  {
4838  ostringstream os;
4839  os << "The sensor is outside (or at the limit) of the model "
4840  << "atmosphere but\nlooks in the wrong direction (wrong sign "
4841  << "for the zenith angle?).\nThis case includes nadir "
4842  << "looking exactly at the latitude end points.";
4843  throw runtime_error( os.str() );
4844  }
4845 
4846  // We can here set ppc and n as we are outside the atmosphere
4847  ppath.nreal = 1.0;
4848  ppath.ngroup = 1.0;
4849  const Numeric r_p = r_e + rte_pos[0];
4850  ppath.constant = geometrical_ppc( r_p, rte_los[0] );
4851 
4852  // Determine TOA radii, and min and max value
4853  Vector r_toa(llat+1);
4854  Numeric r_toa_min = 99e99, r_toa_max = -1;
4855  for( Index ilat=0; ilat<=llat; ilat++ )
4856  {
4857  r_toa[ilat] = refell2r( refellipsoid, lat_grid[ilat] ) +
4858  z_field(lp,ilat,0);
4859  if( r_toa[ilat] < r_toa_min )
4860  { r_toa_min = r_toa[ilat]; }
4861  if( r_toa[ilat] > r_toa_max )
4862  { r_toa_max = r_toa[ilat]; }
4863  }
4864  if( r_p <= r_toa_max )
4865  {
4866  ostringstream os;
4867  os << "The sensor is horizontally outside (or at the limit) of "
4868  << "the model\natmosphere, but is at a radius smaller than "
4869  << "the maximum value of\nthe top-of-the-atmosphere radii. "
4870  << "This is not allowed. Make the\nmodel atmosphere larger "
4871  << "to also cover the sensor position?";
4872  throw runtime_error( os.str() );
4873  }
4874 
4875  // Upward:
4876  if( abs(rte_los[0]) <= 90 )
4877  {
4878  ppath.pos(0,0) = rte_pos[0];
4879  ppath.pos(0,1) = rte_pos[1];
4880  ppath.r[0] = r_e + rte_pos[0];
4881  ppath.los(0,0) = rte_los[0];
4882  //
4883  ppath_set_background( ppath, 1 );
4884  out1 << " ------- WARNING -------: path is totally outside of "
4885  << "the model atmosphere\n";
4886  }
4887 
4888  // Downward:
4889  else
4890  {
4891  bool above=false, ready=false, failed=false;
4892  Numeric rt=-1, latt, lt, lt_old = L_NOT_FOUND;
4893  GridPos gp_latt;
4894  Vector itwt(2);
4895 
4896  // Check if clearly above the model atmosphere
4897  if( ppath.constant >= r_toa_max )
4898  { above=true; ready=true; }
4899  else
4900  { // Otherwise pick out a suitable first test radius
4901  if( islatin || ppath.constant > r_toa_min )
4902  { rt = r_toa_max; }
4903  else
4904  { rt = r_toa_min; }
4905  }
4906 
4907  // Iterate until solution found or moving out from model atm.
4908  //
4909  while( !ready && !failed )
4910  {
4911  // If rt < ppath.constant, ppath above atmosphere
4912  if( rt < ppath.constant )
4913  {
4914  above = true;
4915  ready = true;
4916  }
4917 
4918  else
4919  {
4920  // Calculate length to entrance point at rt
4921  r_crossing_2d( latt, lt, rt, r_p, rte_pos[1], rte_los[0],
4922  ppath.constant );
4923  assert( lt < 9e9 );
4924 
4925  // Entrance outside range of lat_grid = fail
4926  if( latt < lat_grid[0] || latt > lat_grid[llat] )
4927  { failed = true; }
4928 
4929  // OK iteration
4930  else
4931  {
4932  // Converged?
4933  if( abs( lt-lt_old ) < LACC )
4934  { ready = true; }
4935 
4936  // Update rt
4937  lt_old = lt;
4938  gridpos( gp_latt, lat_grid, latt );
4939  interpweights( itwt, gp_latt );
4940  rt = interp( itwt, r_toa, gp_latt );
4941  }
4942  }
4943  } // while
4944 
4945  if( failed )
4946  {
4947  ostringstream os;
4948  os << "The path does not enter the model atmosphere. It "
4949  << "reaches the\ntop of the atmosphere "
4950  << "altitude around latitude " << latt << " deg.";
4951  throw runtime_error( os.str() );
4952  }
4953  else if( above )
4954  {
4955  ppath.pos(0,0) = rte_pos[0];
4956  ppath.pos(0,1) = rte_pos[1];
4957  ppath.r[0] = r_e + rte_pos[0];
4958  ppath.los(0,0) = rte_los[0];
4959  //
4960  ppath_set_background( ppath, 1 );
4961  out1 << " ------- WARNING -------: path is totally outside "
4962  << "of the model atmosphere\n";
4963  }
4964  else
4965  {
4966  ppath.r[0] = rt;
4967  ppath.pos(0,0) = interp( itwt, z_field(lp,joker,0), gp_latt );
4968  // Calculate za first and use to determine lat
4969  ppath.los(0,0) = geompath_za_at_r( ppath.constant,
4970  rte_los[0], rt );
4971  ppath.pos(0,1) = geompath_lat_at_za( rte_los[0], rte_pos[1],
4972  ppath.los(0,0) );
4973  ppath.end_lstep = lt;
4974 
4975  // Here we know the pressure grid position exactly
4976  ppath.gp_p[0].idx = lp-1;
4977  ppath.gp_p[0].fd[0] = 1;
4978  ppath.gp_p[0].fd[1] = 0;
4979 
4980  // Latitude grid position already calculated
4981  gridpos_copy( ppath.gp_lat[0], gp_latt );
4982 
4983  // Hit with cloudbox reaching TOA?
4984  if( cloudbox_on && cloudbox_limits[1] == lp )
4985  {
4986  Numeric fgp = fractional_gp(gp_latt);
4987  if( fgp >= (Numeric)cloudbox_limits[2] &&
4988  fgp <= (Numeric)cloudbox_limits[3] )
4989  { ppath_set_background( ppath, 3 ); }
4990  }
4991  }
4992  } // Downward
4993  } // Outside
4994  } // End 2D
4995 
4996 
4997  //-- 3D ---------------------------------------------------------------------
4998  else
4999  {
5000  // Index of last latitude and longitude
5001  const Index llat = lat_grid.nelem() - 1;
5002  const Index llon = lon_grid.nelem() - 1;
5003 
5004  // Adjust longitude of rte_pos to range used in lon_grid
5005  Numeric lon2use = rte_pos[2];
5006  resolve_lon( lon2use, lon_grid[0], lon_grid[llon] );
5007 
5008  // End position and LOS
5009  ppath.end_pos[0] = rte_pos[0];
5010  ppath.end_pos[1] = rte_pos[1];
5011  ppath.end_pos[2] = lon2use;
5012  ppath.end_los[0] = rte_los[0];
5013  ppath.end_los[1] = rte_los[1];
5014 
5015  // Is sensor inside range of lat_grid and lon_grid?
5016  // If yes, determine TOA altitude at sensor position
5017  GridPos gp_lat, gp_lon;
5018  Vector itw(4);
5019  bool islatlonin = false;
5020  Numeric r_e; // Ellipsoid radius at sensor position
5021  Numeric z_toa = -99e99;
5022 
5023  if( rte_pos[1] >= lat_grid[0] && rte_pos[1] <= lat_grid[llat] &&
5024  lon2use >= lon_grid[0] && lon2use <= lon_grid[llon] )
5025  {
5026  islatlonin = true;
5027  gridpos( gp_lat, lat_grid, rte_pos[1] );
5028  gridpos( gp_lon, lon_grid, lon2use );
5029  interpweights( itw, gp_lat, gp_lon );
5030  z_toa = interp( itw, z_field(lp,joker,joker), gp_lat, gp_lon );
5031  r_e = refell2d( refellipsoid, lat_grid, gp_lat );
5032  }
5033  else
5034  { r_e = refell2r( refellipsoid, rte_pos[1] ); }
5035 
5036  // Sensor is inside the model atmosphere:
5037  if( islatlonin && rte_pos[0] < z_toa )
5038  {
5039  const Numeric z_s = interp( itw, z_surface, gp_lat, gp_lon );
5040 
5041  // Check that the sensor is above the surface
5042  if( (rte_pos[0] + RTOL) < z_s )
5043  {
5044  ostringstream os;
5045  os << "The ppath starting point is placed "
5046  << (z_s - rte_pos[0])/1e3 << " km below the surface.";
5047  throw runtime_error(os.str());
5048  }
5049 
5050  // Set ppath
5051  ppath.pos(0,joker) = ppath.end_pos;
5052  ppath.r[0] = r_e + rte_pos[0];
5053  ppath.los(0,joker) = ppath.end_los;
5054 
5055  // Determine gp_p (gp_lat and gp_lon are recalculated, but ...)
5056  rte_pos2gridpos( ppath.gp_p[0], ppath.gp_lat[0], ppath.gp_lon[0],
5057  atmosphere_dim, p_grid, lat_grid, lon_grid, z_field, rte_pos );
5058  gridpos_check_fd( ppath.gp_p[0] );
5059  gridpos_check_fd( ppath.gp_lat[0] );
5060  gridpos_check_fd( ppath.gp_lon[0] );
5061 
5062  // Is the sensor on the surface looking down?
5063  // If yes and the sensor is inside the cloudbox, the background will
5064  // be changed below.
5065  if( ppath.pos(0,0) <= z_s )
5066  {
5067  // Calculate radial slope of the surface
5068  Numeric c1, c2;
5069  plevel_slope_3d( c1, c2, lat_grid, lon_grid, refellipsoid,
5070  z_surface, gp_lat, gp_lon, ppath.los(0,1) );
5071 
5072  // Calculate angular tilt of the surface
5073  const Numeric atilt = plevel_angletilt( r_e + z_s, c1 );
5074 
5075  // Are we looking down into the surface?
5076  // If yes and the sensor is inside the cloudbox, the background
5077  // will be changed below.
5078  if( is_los_downwards( ppath.los(0,0), atilt ) )
5079  { ppath_set_background( ppath, 2 ); }
5080  }
5081 
5082  // Outside cloudbox:
5083  // Check sensor position with respect to cloud box.
5084  if( cloudbox_on && !ppath_inside_cloudbox_do )
5085  {
5086  const Numeric fgp = fractional_gp( ppath.gp_p[0] );
5087  const Numeric fgl = fractional_gp( ppath.gp_lat[0] );
5088  const Numeric fgo = fractional_gp( ppath.gp_lon[0] );
5089  if( fgp >= cloudbox_limits[0] && fgp <= cloudbox_limits[1] &&
5090  fgl >= cloudbox_limits[2] && fgl <= cloudbox_limits[3] &&
5091  fgo >= cloudbox_limits[4] && fgo <= cloudbox_limits[5] )
5092  { ppath_set_background( ppath, 4 ); }
5093  }
5094 
5095  // Inside cloudbox:
5096  DEBUG_ONLY( if( ppath_inside_cloudbox_do )
5097  {
5098  const Numeric fgp = fractional_gp( ppath.gp_p[0] );
5099  const Numeric fgl = fractional_gp( ppath.gp_lat[0] );
5100  const Numeric fgo = fractional_gp( ppath.gp_lon[0] );
5101  assert( fgp >= cloudbox_limits[0] && fgp <= cloudbox_limits[1] &&
5102  fgl >= cloudbox_limits[2] && fgl <= cloudbox_limits[3] &&
5103  fgo >= cloudbox_limits[4] && fgo <= cloudbox_limits[5] );
5104  } )
5105  }
5106 
5107  // Sensor is outside the model atmosphere:
5108  else
5109  {
5110  // Handle cases when the sensor appears to look the wrong way in
5111  // the north-south direction
5112  if( ( rte_pos[1] <= lat_grid[0] && abs( rte_los[1] ) >= 90 ) ||
5113  ( rte_pos[1] >= lat_grid[llat] && abs( rte_los[1] ) <= 90 ) )
5114  {
5115  ostringstream os;
5116  os << "The sensor is north or south (or at the limit) of the "
5117  << "model atmosphere\nbut looks in the wrong direction.";
5118  throw runtime_error( os.str() );
5119  }
5120 
5121  // Handle cases when the sensor appears to look the wrong way in
5122  // the west-east direction. We demand that the sensor is inside the
5123  // range of lon_grid even if all longitudes are covered.
5124  if( ( lon2use <= lon_grid[0] && rte_los[1] < 0 ) ||
5125  ( lon2use >= lon_grid[llon] && rte_los[1] > 0 ) )
5126  {
5127  ostringstream os;
5128  os << "The sensor is east or west (or at the limit) of the "
5129  << "model atmosphere\nbut looks in the wrong direction.";
5130  throw runtime_error( os.str() );
5131  }
5132 
5133  // We can here set ppc and n as we are outside the atmosphere
5134  ppath.nreal = 1.0;
5135  ppath.ngroup = 1.0;
5136  const Numeric r_p = r_e + rte_pos[0];
5137  ppath.constant = geometrical_ppc( r_p, rte_los[0] );
5138 
5139  // Determine TOA radii, and min and max value
5140  Matrix r_toa(llat+1,llon+1);
5141  Numeric r_toa_min = 99e99, r_toa_max = -1;
5142  for( Index ilat=0; ilat<=llat; ilat++ )
5143  {
5144  const Numeric r_lat = refell2r(refellipsoid,lat_grid[ilat]);
5145  for( Index ilon=0; ilon<=llon; ilon++ )
5146  {
5147  r_toa(ilat,ilon) = r_lat+ z_field(lp,ilat,ilon);
5148  if( r_toa(ilat,ilon) < r_toa_min )
5149  { r_toa_min = r_toa(ilat,ilon); }
5150  if( r_toa(ilat,ilon) > r_toa_max )
5151  { r_toa_max = r_toa(ilat,ilon); }
5152  }
5153  }
5154 
5155  if( r_p <= r_toa_max )
5156  {
5157  ostringstream os;
5158  os << "The sensor is horizontally outside (or at the limit) of "
5159  << "the model\natmosphere, but is at a radius smaller than "
5160  << "the maximum value of\nthe top-of-the-atmosphere radii. "
5161  << "This is not allowed. Make the\nmodel atmosphere larger "
5162  << "to also cover the sensor position?";
5163  throw runtime_error( os.str() );
5164  }
5165 
5166  // Upward:
5167  if( rte_los[0] <= 90 )
5168  {
5169  ppath.pos(0,0) = rte_pos[0];
5170  ppath.pos(0,1) = rte_pos[1];
5171  ppath.pos(0,1) = lon2use;
5172  ppath.r[0] = r_e + rte_pos[0];
5173  ppath.los(0,0) = rte_los[0];
5174  ppath.los(0,1) = rte_los[1];
5175  //
5176  ppath_set_background( ppath, 1 );
5177  out1 << " ------- WARNING -------: path is totally outside of "
5178  << "the model atmosphere\n";
5179  }
5180 
5181  // Downward:
5182  else
5183  {
5184  bool above=false, ready=false, failed=false;
5185  Numeric rt=-1, latt, lont, lt, lt_old = L_NOT_FOUND;
5186  GridPos gp_latt, gp_lont;
5187  Vector itwt(4);
5188 
5189  // Check if clearly above the model atmosphere
5190  if( ppath.constant >= r_toa_max )
5191  { above=true; ready=true; }
5192  else
5193  { // Otherwise pick out a suitable first test radius
5194  if( islatlonin || ppath.constant > r_toa_min )
5195  { rt = r_toa_max; }
5196  else
5197  { rt = r_toa_min; }
5198  }
5199 
5200  // Sensor pos and LOS in cartesian coordinates
5201  Numeric x, y, z, dx, dy, dz;
5202  poslos2cart( x, y, z, dx, dy, dz, r_p, rte_pos[1], lon2use,
5203  rte_los[0], rte_los[1] );
5204 
5205  // Iterate until solution found or moving out from model atm.
5206  //
5207  while( !ready && !failed )
5208  {
5209  // If rt < ppath.constant, ppath above atmosphere
5210  if( rt < ppath.constant )
5211  {
5212  above = true;
5213  ready = true;
5214  }
5215 
5216  else
5217  {
5218  // Calculate lat and lon for entrance point at rt
5219  r_crossing_3d( latt, lont, lt, rt, r_p, rte_pos[1],
5220  lon2use, rte_los[0], ppath.constant,
5221  x, y, z, dx, dy, dz );
5222  resolve_lon( lont, lon_grid[0], lon_grid[llon] );
5223 
5224  // Entrance outside range of lat/lon_grids = fail
5225  if( latt < lat_grid[0] || latt > lat_grid[llat] ||
5226  lont < lon_grid[0] || lont > lon_grid[llon] )
5227  { failed = true; }
5228 
5229  // OK iteration
5230  else
5231  {
5232  // Converged?
5233  if( abs( lt-lt_old ) < LACC )
5234  { ready = true; }
5235 
5236  // Update rt
5237  lt_old = lt;
5238  gridpos( gp_latt, lat_grid, latt );
5239  gridpos( gp_lont, lon_grid, lont );
5240  interpweights( itwt, gp_latt, gp_lont );
5241  rt = interp( itwt, r_toa, gp_latt, gp_lont );
5242  }
5243  }
5244  } // while
5245 
5246  if( failed )
5247  {
5248  ostringstream os;
5249  os << "The path does not enter the model atmosphere. It\n"
5250  << "reaches the top of the atmosphere altitude around:\n"
5251  << " lat: " << latt << " deg.\n lon: " << lont
5252  << " deg.";
5253  throw runtime_error( os.str() );
5254  }
5255  else if( above )
5256  {
5257  ppath.pos(0,0) = rte_pos[0];
5258  ppath.pos(0,1) = rte_pos[1];
5259  ppath.pos(0,1) = lon2use;
5260  ppath.r[0] = r_e + rte_pos[0];
5261  ppath.los(0,0) = rte_los[0];
5262  ppath.los(0,1) = rte_los[1];
5263  //
5264  ppath_set_background( ppath, 1 );
5265  out1 << " ------- WARNING -------: path is totally outside "
5266  << "of the model atmosphere\n";
5267  }
5268  else
5269  {
5270  // Calculate lt for last rt, and use it to determine pos/los
5271  lt = geompath_l_at_r( ppath.constant, r_p ) -
5272  geompath_l_at_r( ppath.constant, rt );
5273  cart2poslos( ppath.r[0], ppath.pos(0,1), ppath.pos(0,2),
5274  ppath.los(0,0), ppath.los(0,1),x+dx*lt, y+dy*lt,
5275  z+dz*lt, dx, dy, dz, ppath.constant, rte_pos[1],
5276  lon2use, rte_los[0], rte_los[1] );
5277  assert( abs( ppath.r[0] -rt ) < RTOL );
5278  resolve_lon( ppath.pos(0,2), lon_grid[0], lon_grid[llon] );
5279  //
5280  ppath.pos(0,0) = interp( itwt, z_field(lp,joker,joker),
5281  gp_latt, gp_lont );
5282  ppath.end_lstep = lt;
5283 
5284  // Here we know the pressure grid position exactly
5285  ppath.gp_p[0].idx = lp-1;
5286  ppath.gp_p[0].fd[0] = 1;
5287  ppath.gp_p[0].fd[1] = 0;
5288 
5289  // Lat and lon grid position already calculated
5290  gridpos_copy( ppath.gp_lat[0], gp_latt );
5291  gridpos_copy( ppath.gp_lon[0], gp_lont );
5292 
5293  // Hit with cloudbox reaching TOA?
5294  if( cloudbox_on && cloudbox_limits[1] == lp )
5295  {
5296  Numeric fgp1 = fractional_gp(gp_latt);
5297  Numeric fgp2 = fractional_gp(gp_lont);
5298  if( fgp1 >= (Numeric)cloudbox_limits[2] &&
5299  fgp1 <= (Numeric)cloudbox_limits[3] &&
5300  fgp2 >= (Numeric)cloudbox_limits[4] &&
5301  fgp2 <= (Numeric)cloudbox_limits[5])
5302  { ppath_set_background( ppath, 3 ); }
5303  }
5304  }
5305  } // Downward
5306  } // Outside
5307  } // End 3D
5308 }
5309 
5310 
5311 
5312 
5313 
5315 
5345  Workspace& ws,
5346  Ppath& ppath,
5347  const Agenda& ppath_step_agenda,
5348  const Index& atmosphere_dim,
5349  const Vector& p_grid,
5350  const Vector& lat_grid,
5351  const Vector& lon_grid,
5352  const Tensor3& t_field,
5353  const Tensor3& z_field,
5354  const Tensor4& vmr_field,
5355  const Vector& f_grid,
5356  const Vector& refellipsoid,
5357  const Matrix& z_surface,
5358  const Index& cloudbox_on,
5359  const ArrayOfIndex& cloudbox_limits,
5360  const Vector& rte_pos,
5361  const Vector& rte_los,
5362  const Numeric& ppath_lraytrace,
5363  const bool& ppath_inside_cloudbox_do,
5364  const Verbosity& verbosity)
5365 {
5366  // This function is a WSM but it is normally only called from yCalc.
5367  // For that reason, this function does not repeat input checks that are
5368  // performed in yCalc, it only performs checks regarding the sensor
5369  // position and LOS.
5370 
5371  //--- Check input -----------------------------------------------------------
5372  chk_rte_pos( atmosphere_dim, rte_pos );
5373  chk_rte_los( atmosphere_dim, rte_los );
5374  if( ppath_inside_cloudbox_do && !cloudbox_on )
5375  throw runtime_error( "The WSV *ppath_inside_cloudbox_do* can only be set "
5376  "to 1 if also *cloudbox_on* is 1." );
5377  //--- End: Check input ------------------------------------------------------
5378 
5379 
5380  // Initiate the partial Ppath structure.
5381  // The function doing the work sets ppath_step to the point of the path
5382  // inside the atmosphere closest to the sensor, if the path is at all inside
5383  // the atmosphere.
5384  // If the background field is set by the function this flags that there is no
5385  // path to follow (for example when the sensor is inside the cloud box).
5386  // The function checks also that the sensor position and LOS give an
5387  // allowed path.
5388  //
5389  Ppath ppath_step;
5390  //
5391  ppath_start_stepping( ppath_step, atmosphere_dim, p_grid, lat_grid,
5392  lon_grid, z_field, refellipsoid, z_surface,
5393  cloudbox_on, cloudbox_limits, ppath_inside_cloudbox_do,
5394  rte_pos, rte_los, verbosity );
5395  // For debugging:
5396  //Print( ppath_step, 0, verbosity );
5397 
5398  // The only data we need to store from this initial ppath_step is
5399  // end_pos/los/lstep
5400  const Numeric end_lstep = ppath_step.end_lstep;
5401  const Vector end_pos = ppath_step.end_pos;
5402  const Vector end_los = ppath_step.end_los;
5403 
5404  // Perform propagation path steps until the starting point is found, which
5405  // is flagged by ppath_step by setting the background field.
5406  //
5407  // The results of each step, returned by ppath_step_agenda as a new
5408  // ppath_step, are stored as an array of Ppath structures.
5409  //
5410  Array<Ppath> ppath_array(0);
5411  Index np = 1; // Counter for number of points of the path
5412  Index istep = 0; // Counter for number of steps
5413  //
5414  const Index imax_p = p_grid.nelem() - 1;
5415  const Index imax_lat = lat_grid.nelem() - 1;
5416  const Index imax_lon = lon_grid.nelem() - 1;
5417  //
5418  bool ready = ppath_what_background( ppath_step );
5419  //
5420  while( !ready )
5421  {
5422  // Call ppath_step agenda.
5423  // The new path step is added to *ppath_array* last in the while block
5424  //
5425  istep++;
5426  //
5427  ppath_step_agendaExecute( ws, ppath_step, ppath_lraytrace, t_field,
5428  z_field, vmr_field, f_grid,
5429  ppath_step_agenda );
5430  // For debugging:
5431  //Print( ppath_step, 0, verbosity );
5432 
5433  // Number of points in returned path step
5434  const Index n = ppath_step.np;
5435 
5436  // Increase the total number
5437  np += n - 1;
5438 
5439  if( istep > 1e4 )
5440  throw runtime_error(
5441  "10 000 path points have been reached. Is this an infinite loop?" );
5442 
5443 
5444  //----------------------------------------------------------------------
5445  //--- Check if some boundary is reached
5446  //----------------------------------------------------------------------
5447 
5448  //--- Outside cloud box ------------------------------------------------
5449  if( !ppath_inside_cloudbox_do )
5450  {
5451  // Check if the top of the atmosphere is reached
5452  // (Perform first a strict check, and if upward propagation also a
5453  // non-strict check. The later to avoid fatal "fixes" at TOA.)
5454  if( is_gridpos_at_index_i( ppath_step.gp_p[n-1], imax_p ) ||
5455  ( abs(ppath_step.los(n-1,0))<90 &&
5456  is_gridpos_at_index_i( ppath_step.gp_p[n-1], imax_p, false )))
5457  {
5458  ppath_set_background( ppath_step, 1 );
5459  }
5460 
5461  // Check that path does not exit at a latitude or longitude end face
5462  if( atmosphere_dim == 2 )
5463  {
5464  // Latitude
5465  if( is_gridpos_at_index_i( ppath_step.gp_lat[n-1], 0 ) )
5466  {
5467  ostringstream os;
5468  os << "The path exits the atmosphere through the lower "
5469  << "latitude end face.\nThe exit point is at an altitude"
5470  << " of " << ppath_step.pos(n-1,0)/1e3 << " km.";
5471  throw runtime_error( os.str() );
5472  }
5473  if( is_gridpos_at_index_i( ppath_step.gp_lat[n-1], imax_lat ) )
5474  {
5475  ostringstream os;
5476  os << "The path exits the atmosphere through the upper "
5477  << "latitude end face.\nThe exit point is at an altitude"
5478  << " of " << ppath_step.pos(n-1,0)/1e3 << " km.";
5479  throw runtime_error( os.str() );
5480  }
5481  }
5482  if( atmosphere_dim == 3 )
5483  {
5484  // Latitude
5485  if( lat_grid[0] > -90 &&
5486  is_gridpos_at_index_i( ppath_step.gp_lat[n-1], 0 ) )
5487  {
5488  ostringstream os;
5489  os << "The path exits the atmosphere through the lower "
5490  << "latitude end face.\nThe exit point is at an altitude"
5491  << " of " << ppath_step.pos(n-1,0)/1e3 << " km.";
5492  throw runtime_error( os.str() );
5493  }
5494  if( lat_grid[imax_lat] < 90 &&
5495  is_gridpos_at_index_i( ppath_step.gp_lat[n-1], imax_lat ) )
5496  {
5497  ostringstream os;
5498  os << "The path exits the atmosphere through the upper "
5499  << "latitude end face.\nThe exit point is at an altitude"
5500  << " of " << ppath_step.pos(n-1,0)/1e3 << " km.";
5501  throw runtime_error( os.str() );
5502  }
5503 
5504  // Longitude
5505  // Note that it must be if and else if here. Otherwise e.g. -180
5506  // will be shifted to 180 and then later back to -180.
5507  if( is_gridpos_at_index_i( ppath_step.gp_lon[n-1], 0 ) &&
5508  ppath_step.los(n-1,1) < 0 &&
5509  abs( ppath_step.pos(n-1,1) ) < 90 )
5510  {
5511  // Check if the longitude point can be shifted +360 degrees
5512  if( lon_grid[imax_lon] - lon_grid[0] >= 360 )
5513  {
5514  ppath_step.pos(n-1,2) = ppath_step.pos(n-1,2) + 360;
5515  gridpos( ppath_step.gp_lon[n-1], lon_grid,
5516  ppath_step.pos(n-1,2) );
5517  }
5518  else
5519  {
5520  ostringstream os;
5521  os << "The path exits the atmosphere through the lower "
5522  << "longitude end face.\nThe exit point is at an "
5523  << "altitude of " << ppath_step.pos(n-1,0)/1e3
5524  << " km.";
5525  throw runtime_error( os.str() );
5526  }
5527  }
5528  else if(
5529  is_gridpos_at_index_i( ppath_step.gp_lon[n-1], imax_lon ) &&
5530  ppath_step.los(n-1,1) > 0 &&
5531  abs( ppath_step.pos(n-1,1) ) < 90 )
5532  {
5533  // Check if the longitude point can be shifted -360 degrees
5534  if( lon_grid[imax_lon] - lon_grid[0] >= 360 )
5535  {
5536  ppath_step.pos(n-1,2) = ppath_step.pos(n-1,2) - 360;
5537  gridpos( ppath_step.gp_lon[n-1], lon_grid,
5538  ppath_step.pos(n-1,2) );
5539  }
5540  else
5541  {
5542  ostringstream os;
5543  os << "The path exits the atmosphere through the upper "
5544  << "longitude end face.\nThe exit point is at an "
5545  << "altitude of " << ppath_step.pos(n-1,0)/1e3
5546  << " km.";
5547  throw runtime_error( os.str() );
5548  }
5549  }
5550  }
5551 
5552 
5553  // Check if there is an intersection with an active cloud box
5554  if( cloudbox_on )
5555  {
5556  Numeric ipos = fractional_gp( ppath_step.gp_p[n-1] );
5557 
5558  if( ipos >= Numeric( cloudbox_limits[0] ) &&
5559  ipos <= Numeric( cloudbox_limits[1] ) )
5560  {
5561  if( atmosphere_dim == 1 )
5562  { ppath_set_background( ppath_step, 3 ); }
5563  else
5564  {
5565  ipos = fractional_gp( ppath_step.gp_lat[n-1] );
5566 
5567  if( ipos >= Numeric( cloudbox_limits[2] ) &&
5568  ipos <= Numeric( cloudbox_limits[3] ) )
5569  {
5570  if( atmosphere_dim == 2 )
5571  { ppath_set_background( ppath_step, 3 ); }
5572  else
5573  {
5574 
5575  ipos = fractional_gp( ppath_step.gp_lon[n-1] );
5576 
5577  if( ipos >= Numeric( cloudbox_limits[4] ) &&
5578  ipos <= Numeric( cloudbox_limits[5] ) )
5579  { ppath_set_background( ppath_step, 3 ); }
5580  }
5581  }
5582  }
5583  }
5584  }
5585  }
5586 
5587  //--- Inside cloud box -------------------------------------------------
5588  else
5589  {
5590  // A first version just checked if point was at or outside any
5591  // boundary but numerical problems could cause that the start point
5592  // was taken as the exit point. So check of ppath direction had to be
5593  // added. Fractional distances used for this.
5594 
5595  // Pressure dimension
5596  Numeric ipos1 = fractional_gp( ppath_step.gp_p[n-1] );
5597  Numeric ipos2 = fractional_gp( ppath_step.gp_p[n-2] );
5598  assert( ipos1 >= cloudbox_limits[0] );
5599  assert( ipos1 <= cloudbox_limits[1] );
5600  if( ipos1 <= cloudbox_limits[0] && ipos1 < ipos2 )
5601  { ppath_set_background( ppath_step, 3 ); }
5602 
5603  else if( ipos1 >= Numeric( cloudbox_limits[1] ) && ipos1 > ipos2 )
5604  { ppath_set_background( ppath_step, 3 ); }
5605 
5606  else if( atmosphere_dim > 1 )
5607  {
5608  // Latitude dimension
5609  ipos1 = fractional_gp( ppath_step.gp_lat[n-1] );
5610  ipos2 = fractional_gp( ppath_step.gp_lat[n-2] );
5611  assert( ipos1 >= cloudbox_limits[2] );
5612  assert( ipos1 <= cloudbox_limits[3] );
5613  if( ipos1 <= Numeric( cloudbox_limits[2] ) && ipos1 < ipos2 )
5614  { ppath_set_background( ppath_step, 3 ); }
5615 
5616  else if( ipos1 >= Numeric( cloudbox_limits[3] ) && ipos1 > ipos2 )
5617  { ppath_set_background( ppath_step, 3 ); }
5618 
5619  else if ( atmosphere_dim > 2 )
5620  {
5621  // Longitude dimension
5622  ipos1 = fractional_gp( ppath_step.gp_lon[n-1] );
5623  ipos2 = fractional_gp( ppath_step.gp_lon[n-2] );
5624  assert( ipos1 >= cloudbox_limits[4] );
5625  assert( ipos1 <= cloudbox_limits[5] );
5626  if( ipos1 <= Numeric( cloudbox_limits[4] ) && ipos1<ipos2 )
5627  { ppath_set_background( ppath_step, 3 ); }
5628 
5629  else if( ipos1 >= Numeric( cloudbox_limits[5] ) &&
5630  ipos1 > ipos2 )
5631  { ppath_set_background( ppath_step, 3 ); }
5632  }
5633  }
5634  }
5635  //----------------------------------------------------------------------
5636  //--- End of boundary check
5637  //----------------------------------------------------------------------
5638 
5639  // Ready?
5640  if( ppath_what_background( ppath_step ) )
5641  {
5642  ppath_step.start_pos = ppath_step.pos(n-1,joker);
5643  ppath_step.start_los = ppath_step.los(n-1,joker);
5644  ready = true;
5645  }
5646 
5647  // Put new ppath_step in ppath_array
5648  ppath_array.push_back( ppath_step );
5649 
5650  } // End path steps
5651 
5652  // Combine all structures in ppath_array to form the return Ppath structure.
5653  //
5654  ppath_init_structure( ppath, atmosphere_dim, np );
5655  //
5656  const Index na = ppath_array.nelem();
5657  //
5658  if( na == 0 ) // No path, just the starting point
5659  {
5660  ppath_copy( ppath, ppath_step, 1 );
5661  // To set n for positions inside the atmosphere, ppath_step_agenda
5662  // must be called once. The later is not always the case. A fix to handle
5663  // those cases:
5664  if( ppath_what_background(ppath_step) > 1 )
5665  {
5666  //Print( ppath_step, 0, verbosity );
5667  ppath_step_agendaExecute( ws, ppath_step, ppath_lraytrace, t_field,
5668  z_field, vmr_field, f_grid,
5669  ppath_step_agenda );
5670  ppath.nreal[0] = ppath_step.nreal[0];
5671  ppath.ngroup[0] = ppath_step.ngroup[0];
5672  }
5673  }
5674 
5675  else // Otherwise, merge the array elelments
5676  {
5677  np = 0; // Now used as counter for points moved to ppath
5678  //
5679  for( Index i=0; i<na; i++ )
5680  {
5681  // For the first structure, the first point shall be included.
5682  // For later structures, the first point shall not be included, but
5683  // there will always be at least two points.
5684 
5685  Index n = ppath_array[i].np;
5686 
5687  // First index to include
5688  Index i1 = 1;
5689  if( i == 0 )
5690  { i1 = 0; }
5691 
5692  // Vectors and matrices that can be handled by ranges.
5693  ppath.r[ Range(np,n-i1) ] = ppath_array[i].r[ Range(i1,n-i1) ];
5694  ppath.pos( Range(np,n-i1), joker ) =
5695  ppath_array[i].pos( Range(i1,n-i1), joker );
5696  ppath.los( Range(np,n-i1), joker ) =
5697  ppath_array[i].los( Range(i1,n-i1), joker );
5698  ppath.nreal[ Range(np,n-i1) ] =
5699  ppath_array[i].nreal[ Range(i1,n-i1) ];
5700  ppath.ngroup[ Range(np,n-i1) ] =
5701  ppath_array[i].ngroup[ Range(i1,n-i1) ];
5702  ppath.lstep[ Range(np-i1,n-1) ] = ppath_array[i].lstep;
5703 
5704  // Grid positions must be handled by a loop
5705  for( Index j=i1; j<n; j++ )
5706  { ppath.gp_p[np+j-i1] = ppath_array[i].gp_p[j]; }
5707  if( atmosphere_dim >= 2 )
5708  {
5709  for( Index j=i1; j<n; j++ )
5710  { ppath.gp_lat[np+j-i1] = ppath_array[i].gp_lat[j]; }
5711  }
5712  if( atmosphere_dim == 3 )
5713  {
5714  for( Index j=i1; j<n; j++ )
5715  { ppath.gp_lon[np+j-i1] = ppath_array[i].gp_lon[j]; }
5716  }
5717 
5718  // Increase number of points done
5719  np += n - i1;
5720  }
5721 
5722  // Field just included once:
5723  // end_pos/los/lspace from first path_step (extracted above):
5724  ppath.end_lstep = end_lstep;
5725  ppath.end_pos = end_pos;
5726  ppath.end_los = end_los;
5727  // Constant, background and start_pos/los from last path_step:
5728  ppath.constant = ppath_step.constant;
5729  ppath.background = ppath_step.background;
5730  ppath.start_pos = ppath_step.start_pos;
5731  ppath.start_los = ppath_step.start_los;
5732  ppath.start_lstep = ppath_step.start_lstep;
5733  }
5734 }
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
void latlon_at_aa(Numeric &lat2, Numeric &lon2, const Numeric &lat1, const Numeric &lon1, const Numeric &aa, const Numeric &ddeg)
latlon_at_aa
Definition: geodetic.cc:927
void get_refr_index_2d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat)
get_refr_index_2d
Definition: refraction.cc:221
ArrayOfGridPos gp_lat
Definition: ppath.h:77
bool is_lon_cyclic(ConstVectorView grid, const Numeric &epsilon)
Check if the given longitude grid is cyclic.
Definition: logic.cc:466
void get_refr_index_3d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat, const Numeric &lon)
Definition: refraction.cc:315
Numeric geompath_za_at_r(const Numeric &ppc, const Numeric &a_za, const Numeric &r)
geompath_za_at_r
Definition: ppath.cc:147
const Numeric LON_NOT_FOUND
Definition: ppath.cc:95
Numeric constant
Definition: ppath.h:62
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const
Number of elements.
Definition: array.h:176
Declarations having to do with the four output streams.
#define q
Definition: continua.cc:21469
void zaaa2cart(Numeric &dx, Numeric &dy, Numeric &dz, const Numeric &za, const Numeric &aa)
zaaa2cart
Definition: ppath.cc:484
Contains the code to determine roots of polynomials.
Numeric refraction_ppc(const Numeric &r, const Numeric &za, const Numeric &refr_index_air)
refraction_ppc
Definition: ppath.cc:639
void ppath_copy(Ppath &ppath1, const Ppath &ppath2, const Index &ncopy)
ppath_copy
Definition: ppath.cc:1876
Matrix los
Definition: ppath.h:68
void ppath_start_stepping(Ppath &ppath, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const bool &ppath_inside_cloudbox_do, ConstVectorView rte_pos, ConstVectorView rte_los, const Verbosity &verbosity)
ppath_start_stepping
Definition: ppath.cc:4606
void raytrace_3d_linear_basic(Workspace &ws, Array< Numeric > &r_array, Array< Numeric > &lat_array, Array< Numeric > &lon_array, Array< Numeric > &za_array, Array< Numeric > &aa_array, Array< Numeric > &l_array, Array< Numeric > &n_array, Array< Numeric > &ng_array, Index &endface, ConstVectorView refellipsoid, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &lmax, const Agenda &refr_index_air_agenda, const Numeric &lraytrace, const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &rsurface15, const Numeric &rsurface35, const Numeric &rsurface36, const Numeric &rsurface16, const Numeric &r15a, const Numeric &r35a, const Numeric &r36a, const Numeric &r16a, const Numeric &r15b, const Numeric &r35b, const Numeric &r36b, const Numeric &r16b, Numeric r, Numeric lat, Numeric lon, Numeric za, Numeric aa)
raytrace_3d_linear_basic
Definition: ppath.cc:4258
The Vector class.
Definition: matpackI.h:556
Numeric geompath_r_at_l(const Numeric &ppc, const Numeric &l)
geompath_r_at_l
Definition: ppath.cc:275
void do_gridcell_2d_byltest(Vector &r_v, Vector &lat_v, Vector &za_v, Numeric &lstep, Index &endface, const Numeric &r_start0, const Numeric &lat_start0, const Numeric &za_start, const Numeric &l_start, const Index &icall, const Numeric &ppc, const Numeric &lmax, const Numeric &lat1, const Numeric &lat3, const Numeric &r1a, const Numeric &r3a, const Numeric &r3b, const Numeric &r1b, const Numeric &rsurface1, const Numeric &rsurface3)
do_gridcell_2d_byltest
Definition: ppath.cc:2837
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Vector end_pos
Definition: ppath.h:71
void raytrace_2d_linear_basic(Workspace &ws, Array< Numeric > &r_array, Array< Numeric > &lat_array, Array< Numeric > &za_array, Array< Numeric > &l_array, Array< Numeric > &n_array, Array< Numeric > &ng_array, Index &endface, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &lmax, const Agenda &refr_index_air_agenda, const Numeric &lraytrace, const Numeric &lat1, const Numeric &lat3, const Numeric &rsurface1, const Numeric &rsurface3, const Numeric &r1a, const Numeric &r3a, const Numeric &r3b, const Numeric &r1b, Numeric r, Numeric lat, Numeric za)
raytrace_2d_linear_basic
Definition: ppath.cc:3947
Index dim
Definition: ppath.h:60
const Numeric DEG2RAD
void ppath_step_refr_2d(Workspace &ws, Ppath &ppath, ConstVectorView p_grid, ConstVectorView lat_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstVectorView z_surface, const Numeric &lmax, const Agenda &refr_index_air_agenda, const String &rtrace_method, const Numeric &lraytrace)
ppath_step_refr_2d
Definition: ppath.cc:4121
void ppath_step_agendaExecute(Workspace &ws, Ppath &ppath_step, const Numeric ppath_lraytrace, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &f_grid, const Agenda &input_agenda)
Definition: auto_md.cc:14720
The Tensor4 class.
Definition: matpackIV.h:383
The range class.
Definition: matpackI.h:148
Vector lstep
Definition: ppath.h:70
void rotationmat3D(Matrix &R, ConstVectorView vrot, const Numeric &a)
rotationmat3D
Definition: ppath.cc:520
void ppath_step_geom_3d(Ppath &ppath, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Numeric &lmax)
ppath_step_geom_3d
Definition: ppath.cc:3553
void cart2zaaa(Numeric &za, Numeric &aa, const Numeric &dx, const Numeric &dy, const Numeric &dz)
cart2zaaa
Definition: ppath.cc:443
void ppath_end_2d(Ppath &ppath, ConstVectorView r_v, ConstVectorView lat_v, ConstVectorView za_v, ConstVectorView lstep, ConstVectorView n_v, ConstVectorView ng_v, ConstVectorView lat_grid, ConstMatrixView z_field, ConstVectorView refellipsoid, const Index &ip, const Index &ilat, const Index &endface, const Numeric &ppc)
ppath_end_2d
Definition: ppath.cc:2240
Numeric fractional_gp(const GridPos &gp)
fractional_gp
const Numeric RAD2DEG
void chk_rte_los(const Index &atmosphere_dim, ConstVectorView rte_los)
chk_rte_los
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:679
Matrix pos
Definition: ppath.h:67
void plevel_slope_3d(Numeric &c1, Numeric &c2, const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &r15, const Numeric &r35, const Numeric &r36, const Numeric &r16, const Numeric &lat, const Numeric &lon, const Numeric &aa)
plevel_slope_3d
Definition: ppath.cc:1354
Vector ngroup
Definition: ppath.h:75
void raytrace_1d_linear_basic(Workspace &ws, Array< Numeric > &r_array, Array< Numeric > &lat_array, Array< Numeric > &za_array, Array< Numeric > &l_array, Array< Numeric > &n_array, Array< Numeric > &ng_array, Index &endface, ConstVectorView p_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &lmax, const Agenda &refr_index_air_agenda, const Numeric &lraytrace, const Numeric &rsurface, const Numeric &r1, const Numeric &r3, Numeric r, Numeric lat, Numeric za)
raytrace_1d_linear_basic
Definition: ppath.cc:3658
Numeric refell2d(ConstVectorView refellipsoid, ConstVectorView lat_grid, const GridPos gp)
refell2d
Definition: geodetic.cc:1244
Index gridpos2gridrange(const GridPos &gp, const bool &upwards)
gridpos2gridrange
Vector start_pos
Definition: ppath.h:64
Vector r
Definition: ppath.h:69
Numeric rslope_crossing3d(const Numeric &rp, const Numeric &za, const Numeric &r0, Numeric c1, Numeric c2)
rslope_crossing3d
Definition: ppath.cc:1511
const Numeric LAT_NOT_FOUND
Definition: ppath.cc:94
void get_refr_index_1d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r)
get_refr_index_1d
Definition: refraction.cc:141
Numeric geompath_lat_at_za(const Numeric &za0, const Numeric &lat0, const Numeric &za)
geompath_lat_at_za
Definition: ppath.cc:218
bool is_gridpos_at_index_i(const GridPos &gp, const Index &i, const bool &strict)
is_gridpos_at_index_i
void resolve_lon(Numeric &lon, const Numeric &lon5, const Numeric &lon6)
resolve_lon
Definition: ppath.cc:675
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Numeric geometrical_ppc(const Numeric &r, const Numeric &za)
geometrical_ppc
Definition: ppath.cc:118
Structure to store a grid position.
Definition: interpolation.h:74
A constant view of a Tensor4.
Definition: matpackIV.h:141
Index ppath_what_background(const Ppath &ppath)
ppath_what_background
Definition: ppath.cc:1835
Numeric end_lstep
Definition: ppath.h:73
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:473
const Numeric L_NOT_FOUND
Definition: ppath.cc:93
This file contains the definition of Array.
void gridpos_check_fd(GridPos &gp)
gridpos_check_fd
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
void map_daa(Numeric &za, Numeric &aa, const Numeric &za0, const Numeric &aa0, const Numeric &aa_grid)
Definition: ppath.cc:576
void plevel_slope_2d(Numeric &c1, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstVectorView z_surf, const GridPos &gp, const Numeric &za)
plevel_slope_2d
Definition: ppath.cc:782
Vector end_los
Definition: ppath.h:72
#define dmin(a, b)
Definition: continua.cc:20462
The implementation for String, the ARTS string class.
Definition: mystring.h:63
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
void ppath_set_background(Ppath &ppath, const Index &case_nr)
ppath_set_background
Definition: ppath.cc:1790
The Tensor3 class.
Definition: matpackIII.h:348
String background
Definition: ppath.h:63
Numeric rslope_crossing2d(const Numeric &rp, const Numeric &za, const Numeric &r0, Numeric c1)
rslope_crossing2d
Definition: ppath.cc:992
void do_gridrange_1d(Vector &r_v, Vector &lat_v, Vector &za_v, Numeric &lstep, Index &endface, const Numeric &r_start0, const Numeric &lat_start, const Numeric &za_start, const Numeric &ppc, const Numeric &lmax, const Numeric &ra, const Numeric &rb, const Numeric &rsurface)
do_gridrange_1d
Definition: ppath.cc:2699
#define DEBUG_ONLY(...)
Definition: debug.h:37
void geompath_from_r1_to_r2(Vector &r, Vector &lat, Vector &za, Numeric &lstep, const Numeric &ppc, const Numeric &r1, const Numeric &lat1, const Numeric &za1, const Numeric &r2, const bool &tanpoint, const Numeric &lmax)
geompath_from_r1_to_r2
Definition: ppath.cc:341
void ppath_end_1d(Ppath &ppath, ConstVectorView r_v, ConstVectorView lat_v, ConstVectorView za_v, ConstVectorView lstep, ConstVectorView n_v, ConstVectorView ng_v, ConstVectorView z_field, ConstVectorView refellipsoid, const Index &ip, const Index &endface, const Numeric &ppc)
ppath_end_1d
Definition: ppath.cc:2052
const Numeric ANGTOL
Definition: ppath.h:103
void ppath_end_3d(Ppath &ppath, ConstVectorView r_v, ConstVectorView lat_v, ConstVectorView lon_v, ConstVectorView za_v, ConstVectorView aa_v, ConstVectorView lstep, ConstVectorView n_v, ConstVectorView ng_v, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, const Index &ip, const Index &ilat, const Index &ilon, const Index &endface, const Numeric &ppc)
ppath_end_3d
Definition: ppath.cc:2527
void chk_rte_pos(const Index &atmosphere_dim, ConstVectorView rte_pos, const bool &is_rte_pos2)
chk_rte_pos
void do_gridcell_3d_byltest(Vector &r_v, Vector &lat_v, Vector &lon_v, Vector &za_v, Vector &aa_v, Numeric &lstep, Index &endface, const Numeric &r_start0, const Numeric &lat_start0, const Numeric &lon_start0, const Numeric &za_start, const Numeric &aa_start, const Numeric &l_start, const Index &icall, const Numeric &ppc, const Numeric &lmax, const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &r15a, const Numeric &r35a, const Numeric &r36a, const Numeric &r16a, const Numeric &r15b, const Numeric &r35b, const Numeric &r36b, const Numeric &r16b, const Numeric &rsurface15, const Numeric &rsurface35, const Numeric &rsurface36, const Numeric &rsurface16)
do_gridcell_3d_byltest
Definition: ppath.cc:3191
Declarations for agendas.
void ppath_init_structure(Ppath &ppath, const Index &atmosphere_dim, const Index &np)
ppath_init_structure
Definition: ppath.cc:1726
#define max(a, b)
Definition: continua.cc:20461
void cart2poslos(Numeric &r, Numeric &lat, Numeric &za, const Numeric &x, const Numeric &z, const Numeric &dx, const Numeric &dz, const Numeric &ppc, const Numeric &lat0, const Numeric &za0)
cart2poslos
Definition: geodetic.cc:125
Numeric start_lstep
Definition: ppath.h:66
#define CREATE_OUT1
Definition: messages.h:212
void ppath_append(Ppath &ppath1, const Ppath &ppath2)
ppath_append
Definition: ppath.cc:1948
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
#define abs(x)
Definition: continua.cc:20458
void ppath_start_1d(Numeric &r_start, Numeric &lat_start, Numeric &za_start, Index &ip, const Ppath &ppath)
ppath_start_1d
Definition: ppath.cc:2017
const Joker joker
const Numeric R_NOT_FOUND
Definition: ppath.cc:92
void find_tanpoint(Index &it, const Ppath ppath)
find_tanpoint
Definition: ppath.cc:705
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
#define dx
Definition: continua.cc:21928
The Matrix class.
Definition: matpackI.h:788
Vector nreal
Definition: ppath.h:74
Numeric plevel_angletilt(const Numeric &r, const Numeric &c1)
plevel_angletilt
Definition: ppath.cc:844
const Numeric RTOL
Definition: ppath.cc:78
void refr_gradients_3d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, Numeric &dndlat, Numeric &dndlon, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat, const Numeric &lon)
refr_gradients_3d
Definition: refraction.cc:557
Numeric refell2r(ConstVectorView refellipsoid, const Numeric &lat)
refell2r
Definition: geodetic.cc:1199
void gridpos_copy(GridPos &gp_new, const GridPos &gp_old)
gridpos_copy
bool is_los_downwards(const Numeric &za, const Numeric &tilt)
is_los_downwards
Definition: ppath.cc:875
Header file for special_interp.cc.
Propagation path structure and functions.
void ppath_step_geom_2d(Ppath &ppath, ConstVectorView lat_grid, ConstMatrixView z_field, ConstVectorView refellipsoid, ConstVectorView z_surface, const Numeric &lmax)
ppath_step_geom_2d
Definition: ppath.cc:3136
void ppath_step_geom_1d(Ppath &ppath, ConstVectorView z_field, ConstVectorView refellipsoid, const Numeric &z_surface, const Numeric &lmax)
ppath_step_geom_1d
Definition: ppath.cc:2783
void ppath_start_3d(Numeric &r_start, Numeric &lat_start, Numeric &lon_start, Numeric &za_start, Numeric &aa_start, Index &ip, Index &ilat, Index &ilon, Numeric &lat1, Numeric &lat3, Numeric &lon5, Numeric &lon6, Numeric &r15a, Numeric &r35a, Numeric &r36a, Numeric &r16a, Numeric &r15b, Numeric &r35b, Numeric &r36b, Numeric &r16b, Numeric &rsurface15, Numeric &rsurface35, Numeric &rsurface36, Numeric &rsurface16, Ppath &ppath, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView refellipsoid, ConstMatrixView z_surface)
ppath_start_3d
Definition: ppath.cc:2353
void r_crossing_3d(Numeric &lat, Numeric &lon, Numeric &l, const Numeric &r_hit, const Numeric &r_start, const Numeric &lat_start, const Numeric &lon_start, const Numeric &za_start, const Numeric &ppc, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz)
r_crossing_3d
Definition: ppath.cc:1638
void r_crossing_2d(Numeric &lat, Numeric &l, const Numeric &r_hit, const Numeric &r_start, const Numeric &lat_start, const Numeric &za_start, const Numeric &ppc)
r_crossing_2d
Definition: ppath.cc:914
Header file for logic.cc.
Numeric rsurf_at_latlon(const Numeric &lat1, const Numeric &lat3, const Numeric &lon5, const Numeric &lon6, const Numeric &r15, const Numeric &r35, const Numeric &r36, const Numeric &r16, const Numeric &lat, const Numeric &lon)
rsurf_at_latlon
Definition: ppath.cc:1294
void adjust_los(VectorView los, const Index &atmosphere_dim)
adjust_los
Definition: rte.cc:81
Numeric geompath_r_at_lat(const Numeric &ppc, const Numeric &lat0, const Numeric &za0, const Numeric &lat)
geompath_r_at_lat
Definition: ppath.cc:299
void ppath_step_refr_1d(Workspace &ws, Ppath &ppath, ConstVectorView p_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, const Numeric &z_surface, const Numeric &lmax, const Agenda &refr_index_air_agenda, const String &rtrace_method, const Numeric &lraytrace)
ppath_step_refr_1d
Definition: ppath.cc:3809
void rte_pos2gridpos(GridPos &gp_p, GridPos &gp_lat, GridPos &gp_lon, const Index &atmosphere_dim, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstVectorView rte_pos)
rte_pos2gridpos
void cart2pol(Numeric &r, Numeric &lat, const Numeric &x, const Numeric &z, const Numeric &lat0, const Numeric &za0)
cart2pol
Definition: geodetic.cc:82
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
#define min(a, b)
Definition: continua.cc:20460
A constant view of a Tensor3.
Definition: matpackIII.h:139
A constant view of a Vector.
Definition: matpackI.h:292
Index np
Definition: ppath.h:61
ArrayOfGridPos gp_lon
Definition: ppath.h:78
void ppath_start_2d(Numeric &r_start, Numeric &lat_start, Numeric &za_start, Index &ip, Index &ilat, Numeric &lat1, Numeric &lat3, Numeric &r1a, Numeric &r3a, Numeric &r3b, Numeric &r1b, Numeric &rsurface1, Numeric &rsurface3, Ppath &ppath, ConstVectorView lat_grid, ConstMatrixView z_field, ConstVectorView refellipsoid, ConstVectorView z_surface)
ppath_start_2d
Definition: ppath.cc:2121
void refr_gradients_1d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r)
refr_gradients_1d
Definition: refraction.cc:412
A constant view of a Matrix.
Definition: matpackI.h:596
void refr_gradients_2d(Workspace &ws, Numeric &refr_index_air, Numeric &refr_index_air_group, Numeric &dndr, Numeric &dndlat, const Agenda &refr_index_air_agenda, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView refellipsoid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, const Numeric &r, const Numeric &lat)
refr_gradients_2d
Definition: refraction.cc:477
Vector start_los
Definition: ppath.h:65
Workspace class.
Definition: workspace_ng.h:47
void cart2sph(Numeric &r, Numeric &lat, Numeric &lon, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &lat0, const Numeric &lon0, const Numeric &za0, const Numeric &aa0)
cart2sph
Definition: geodetic.cc:596
Refraction functions.
#define beta
Definition: continua.cc:21194
Numeric geompath_r_at_za(const Numeric &ppc, const Numeric &za)
geompath_r_at_za
Definition: ppath.cc:190
const Numeric LATLONTOL
Definition: ppath.cc:83
Numeric rsurf_at_lat(const Numeric &lat1, const Numeric &lat3, const Numeric &r1, const Numeric &r3, const Numeric &lat)
rsurf_at_lat
Definition: ppath.cc:743
void ppath_calc(Workspace &ws, Ppath &ppath, const Agenda &ppath_step_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Vector &f_grid, const Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Vector &rte_pos, const Vector &rte_los, const Numeric &ppath_lraytrace, const bool &ppath_inside_cloudbox_do, const Verbosity &verbosity)
ppath_calc
Definition: ppath.cc:5344
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:59
Header file for helper functions for OpenMP.
void ppath_step_refr_3d(Workspace &ws, Ppath &ppath, ConstVectorView p_grid, ConstVectorView lat_grid, ConstVectorView lon_grid, ConstTensor3View z_field, ConstTensor3View t_field, ConstTensor4View vmr_field, ConstVectorView f_grid, ConstVectorView refellipsoid, ConstMatrixView z_surface, const Numeric &lmax, const Agenda &refr_index_air_agenda, const String &rtrace_method, const Numeric &lraytrace)
ppath_step_refr_3d
Definition: ppath.cc:4466
void plevel_crossing_2d(Numeric &r, Numeric &lat, Numeric &l, const Numeric &r_start0, const Numeric &lat_start, const Numeric &za_start, const Numeric &ppc, const Numeric &lat1, const Numeric &lat3, const Numeric &r1, const Numeric &r3, const bool &above)
plevel_crossing_2d
Definition: ppath.cc:1121
const Numeric POLELAT
Definition: ppath.h:94
int poly_root_solve(Matrix &roots, Vector &coeffs)
Definition: poly_roots.cc:100
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
const Numeric LACC
Definition: ppath.cc:88
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition: geodetic.cc:387
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
ArrayOfGridPos gp_p
Definition: ppath.h:76
Numeric geompath_l_at_r(const Numeric &ppc, const Numeric &r)
geompath_l_at_r
Definition: ppath.cc:246
void gridpos_force_end_fd(GridPos &gp, const Index &n)
gridpos_force_end_fd
This file contains the definition of String, the ARTS string class.
Declaration of functions in rte.cc.
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1580