ARTS  2.2.66
ppath_NotUsed.cc
Go to the documentation of this file.
1 
27  double& r,
28  double& lat,
29  double& lon,
30  double& za,
31  double& aa,
32  const double& x,
33  const double& y,
34  const double& z,
35  const double& dx,
36  const double& dy,
37  const double& dz )
38 {
39  // Assert that LOS vector is normalised
40  assert( abs( sqrt( dx*dx + dy*dy + dz*dz ) - 1 ) < 1e-6 );
41 
42  // Spherical coordinates
43  cart2sph( r, lat, lon, x, y, z );
44 
45  // Spherical derivatives
46  const double coslat = cos( DEG2RAD * lat );
47  const double sinlat = sin( DEG2RAD * lat );
48  const double coslon = cos( DEG2RAD * lon );
49  const double sinlon = sin( DEG2RAD * lon );
50  const double dr = coslat*coslon*dx + coslat*sinlon*dy + sinlat*dz;
51  const double dlat = -sinlat*coslon/r*dx - sinlat*sinlon/r*dy + coslat/r*dz;
52  const double dlon = -sinlon/coslat/r*dx + coslon/coslat/r*dy;
53 
54  // LOS angles
55  za = RAD2DEG * acos( dr );
56  //
57  if( za < ANGTOL || za > 180-ANGTOL )
58  { aa = 0; }
59 
60  else if( abs( lat ) <= POLELAT )
61  {
62  aa = RAD2DEG * acos( r * dlat / sin( DEG2RAD * za ) );
63 
64  if( isnan( aa ) )
65  {
66  if( dlat >= 0 )
67  { aa = 0; }
68  else
69  { aa = 180; }
70  }
71  else if( dlon < 0 )
72  { aa = -aa; }
73  }
74 
75  // For lat = +- 90 the azimuth angle gives the longitude along which the
76  // LOS goes
77  else
78  { aa = RAD2DEG * atan2( dy, dx ); }
79 }
80 
81 
83 
96 void cart2sph(
97  double& r,
98  double& lat,
99  double& lon,
100  const double& x,
101  const double& y,
102  const double& z )
103 {
104  r = sqrt( x*x + y*y + z*z );
105  lat = RAD2DEG * asin( z / r );
106  lon = RAD2DEG * atan2( y, x );
107 }
108 
109 
110 
112 
129  const Numeric& r1,
130  const Numeric& lat1,
131  const Numeric& r2,
132  const Numeric& lat2 )
133 {
134  if( lat2 == lat1 )
135  {
136  if( r1 <= r2 )
137  { return 0; }
138  else
139  { return 180; }
140  }
141  else
142  {
143  // Absolute value of the latitude difference
144  const Numeric dlat = abs( lat2 - lat1 );
145 
146  // The zenith angle is obtained by a combination of the lawes of sine
147  // and cosine.
148  Numeric za = dlat + RAD2DEG * asin( r1 * sin( DEG2RAD * dlat ) /
149  sqrt( r1*r1 + r2*r2 - 2 * r1 * r2 * cos( DEG2RAD * dlat ) ) );
150 
151  // Consider viewing direction
152  if( lat2 < lat1 )
153  { za = -za; }
154 
155  return za;
156  }
157 }
158 
159 
161 
191  Numeric& r,
192  Numeric& lon,
193  Numeric& l,
194  const Numeric& lat_hit,
195  const Numeric& lat_start,
196  const Numeric& za_start,
197  const Numeric& x,
198  const Numeric& y,
199  const Numeric& z,
200  const Numeric& dx,
201  const Numeric& dy,
202  const Numeric& dz )
203 {
204  assert( lat_start >= -90 );
205  assert( lat_start <= 90 );
206  assert( lat_hit >= -90 );
207  assert( lat_hit <= 90 );
208 
209  // For za=0/180 and lat+-90 there is no solution
210  if( za_start == 0 || za_start == 180 || abs(lat_hit) == 90 )
211  { l = -1; }
212 
213  // The expressions below can not be used for lat=0
214  else if( abs( lat_hit ) < 1e-7 )
215  { l = -z / dz; }
216 
217  else
218  {
219  const Numeric t2 = pow( tan( DEG2RAD * lat_hit ), 2.0 );
220  const Numeric a = t2 * ( dx*dx + dy*dy ) - dz*dz;
221  const Numeric b = 2 * ( t2 * ( x*dx + y*dy ) - z*dz );
222  const Numeric c = t2 * ( x*x + y*y ) - z*z;
223  const Numeric bb = b * b;
224  const Numeric ac4 = 4 * a * c;
225 
226  // Check if a real solution is possible
227  if( ac4 > bb )
228  { l = -1; }
229  else
230  {
231  const Numeric d = -0.5*b/a;
232  const Numeric e = -0.5*sqrt(b*b-4*a*c)/a;
233  Numeric l1 = d + e;
234  Numeric l2 = d - e;
235 
236  // Both lat and -lat can end up as a solution (the sign is lost as
237  // tan(lat) is squared). A correct solution requires that l>=0 and
238  // that z+l*dz has the same sigh as lat. Set l to -1 if this not
239  // fulfilled.
240  const Numeric zsign = sign( lat_hit );
241  if( l1 > 0 && abs(sign(z+dz*l1)-zsign)>0.01 )
242  { l1 = -1;}
243  if( l2 > 0 && abs(sign(z+dz*l2)-zsign)>0.01 )
244  { l2 = -1;}
245 
246  // If both l1 and l2 are > 0, we want theoretically the smallest
247  // value. However, with lat=lat0 the "zero solution" can deviate
248  // slightly from zero due to numerical issues, and it is not just to
249  // pick the smallest positive value. As a solution, don't accept a
250  // final l below 1e-6 if not both l1 and l2 are inside [0,1e-6].
251  const Numeric lmin = min( l1, l2 );
252  const Numeric lmax = max( l1, l2 );
253  if( lmin >= 0 && lmax < 1e-6 )
254  { l = lmax; }
255  else
256  {
257  if( lmin > 1e-6 )
258  { l = lmin; }
259  else if( lmax > 1e-6 )
260  { l = lmax; }
261  else
262  { l = -1; }
263  }
264  }
265  }
266 
267  if( l > 0 )
268  {
269  const Numeric xp = x+dx*l;
270  const Numeric yp = y+dy*l;
271  r = sqrt( pow(xp,2.0) + pow(yp,2.0) + pow(z+dz*l,2.0) );
272  lon = RAD2DEG * atan2( yp, xp );
273  }
274  else
275  { r = R_NOT_FOUND; lon = LON_NOT_FOUND; l = L_NOT_FOUND; }
276 }
277 
278 
279 
281 
312  Numeric& r,
313  Numeric& lat,
314  Numeric& l,
315  const Numeric& lon_hit,
316  const Numeric& lon_start,
317  const Numeric& za_start,
318  const Numeric& aa_start,
319  const Numeric& x,
320  const Numeric& y,
321  const Numeric& z,
322  const Numeric& dx,
323  const Numeric& dy,
324  const Numeric& dz )
325 {
326  if( lon_hit == lon_start || za_start == 0 || za_start == 180 ||
327  aa_start == 0 || abs(aa_start) == 180 )
328  { l = -1; }
329 
330  else
331  {
332  const double tanlon = tan( DEG2RAD * lon_hit );
333  l = ( y - x*tanlon ) / ( dx*tanlon - dy );
334  }
335 
336  if( l <= 0 )
337  { r = R_NOT_FOUND; lat = LAT_NOT_FOUND; l = L_NOT_FOUND; }
338 
339  else
340  {
341  const Numeric zp = z + dz*l;
342  r = sqrt( pow(x+dx*l,2.0) + pow(y+dy*l,2.0) + pow(zp,2.0) );
343  lat = RAD2DEG * asin( zp / r );
344  }
345 }
346 
347 
348 
350 
385  Numeric& r,
386  Numeric& lat,
387  Numeric& lon,
388  Numeric& l,
389  const Numeric& r_start0,
390  const Numeric& lat_start,
391  const Numeric& lon_start,
392  const Numeric& za_start,
393  const Numeric& aa_start,
394  const Numeric& x,
395  const Numeric& y,
396  const Numeric& z,
397  const Numeric& dx,
398  const Numeric& dy,
399  const Numeric& dz,
400  const Numeric& ppc,
401  const Numeric& lat1,
402  const Numeric& lat3,
403  const Numeric& lon5,
404  const Numeric& lon6,
405  const Numeric& r15,
406  const Numeric& r35,
407  const Numeric& r36,
408  const Numeric& r16,
409  const bool& above )
410 {
411  assert( za_start <= 180 );
412  assert( lat_start >=lat1 && lat_start <= lat3 );
413  assert( lon_start >=lon5 && lon_start <= lon6 );
414 
415  // Zenith case
416  if( za_start < ANGTOL )
417  {
418  if( above )
419  { r = R_NOT_FOUND; lat = LAT_NOT_FOUND;
420  l = L_NOT_FOUND; lon = LON_NOT_FOUND; }
421  else
422  {
423  lat = lat_start;
424  lon = lon_start;
425  r = rsurf_at_latlon( lat1, lat3, lon5, lon6, r15, r35, r36, r16,
426  lat, lon );
427  l = max( 1e-9, r - r_start0 ); // Max to ensure a small positive
428  } // step, to handle numerical issues
429  }
430 
431  // Nadir case
432  else if( za_start > 180-ANGTOL )
433  {
434  if( above )
435  {
436  lat = lat_start;
437  lon = lon_start;
438  r = rsurf_at_latlon( lat1, lat3, lon5, lon6, r15, r35, r36, r16,
439  lat, lon );
440  l = max( 1e-9, r_start0 - r ); // As above
441  }
442  else
443  { r = R_NOT_FOUND; lat = LAT_NOT_FOUND;
444  l = L_NOT_FOUND; lon = LON_NOT_FOUND; }
445  }
446 
447  // The general case
448  else
449  {
450  const Numeric rmin = min( r15, min( r35, min( r36, r16 ) ) );
451  const Numeric rmax = max( r15, max( r35, max( r36, r16 ) ) );
452 
453  // The case of negligible slope
454  if( rmax-rmin < RTOL/10 )
455  {
456  // Set r_start, considering impact of numerical problems
457  Numeric r_start = r_start0;
458  r = r15;
459  if( above )
460  { if( r_start < rmax ) { r_start = r = rmax; } }
461  else
462  { if( r_start > rmin ) { r_start = r = rmin; } }
463 
464  r_crossing_3d( lat, lon, l, r, r_start, lat_start, lon_start,
465  za_start, ppc, x, y, z, dx, dy, dz );
466 
467  // Check if inside [lat1,lat3]
468  if( lat > lat3 || lat < lat1 || lon > lon6 || lon < lon5 )
469  { r = R_NOT_FOUND; lat = LAT_NOT_FOUND; lon = LON_NOT_FOUND; }
470  }
471 
472  // With slope
473  else
474  {
475  // Set r_start, considering impact of numerical problems
476  Numeric r_start = r_start0;
477  if( above )
478  { if( r_start < rmin ) { r_start = rmin; } }
479  else
480  { if( r_start > rmax ) { r_start = rmax; } }
481 
482  // Calculate crossing with closest radius
483  if( r_start > rmax )
484  {
485  r = rmax;
486  r_crossing_3d( lat, lon, l, r, r_start, lat_start, lon_start,
487  za_start, ppc, x, y, z, dx, dy, dz );
488  }
489  else if( r_start < rmin )
490  {
491  r = rmin;
492  r_crossing_3d( lat, lon, l, r, r_start, lat_start, lon_start,
493  za_start, ppc, x, y, z, dx, dy, dz );
494  }
495  else
496  { r = r_start; lat = lat_start; lon = lon_start; l = 0; }
497 
498  // lat/lon must be inside if relevant to continue
499  if( lat < lat1 || lat > lat3 || lon < lon5 || lon > lon6 )
500  { r = R_NOT_FOUND; } // lat and lon already set by r_crossing_3d
501 
502  // Otherwise continue from found point, considering the level slope
503  else
504  {
505  // Level radius at lat/lon
506  const Numeric rpl = rsurf_at_latlon( lat1, lat3, lon5, lon6,
507  r15, r35, r36, r16, lat, lon );
508 
509  // Make adjustment if numerical problems
510  if( above )
511  { if( r < rpl ) { r = rpl; } }
512  else
513  { if( r > rpl ) { r = rpl; } }
514 
515  // Azimuth angle:
516  Numeric d1, d2, d3, za, aa;
517  if( l == 0 )
518  { za = za_start; aa = aa_start; }
519  else
520  { cart2poslos( d1, d2, d3, za, aa, x+dx*l, y+dy*l, z+dz*l,
521  dx, dy, dz, ppc, lat_start, lon_start,
522  za_start, aa_start );
523  assert( abs(d1-r) < 1e-3 );
524  assert( abs(d2-lat) < 1e-8 );
525  assert( abs(d3-lon) < 1e-8 );
526  }
527 
528  // Level slope at lat/lon
529  Numeric c1, c2;
530  plevel_slope_3d( c1, c2, lat1, lat3, lon5, lon6,
531  r15, r35, r36, r16, lat, lon, aa );
532 
533  // Angular distance from present point to actual crossing
534  const Numeric dang = rslope_crossing3d( r, za, rpl, c1, c2 );
535 
536  // Lat and lon at dang
537  Numeric lat2, lon2;
538  //
539  latlon_at_aa( lat2, lon2, lat, lon, aa, dang );
540  //
541  lat = lat2;
542  lon = lon2; resolve_lon( lon, lon5, lon6 );
543 
544  // lat/lon still inside gridbox? If yes, update r and l
545  if( lat < lat1 || lat > lat3 || lon < lon5 || lon > lon6 )
546  { r = R_NOT_FOUND; lat = LAT_NOT_FOUND;
547  l = L_NOT_FOUND; lon = LON_NOT_FOUND; }
548  else
549  {
550  r = rsurf_at_latlon( lat1, lat3, lon5, lon6,
551  r15, r35, r36, r16, lat, lon );
552  distance3D( l, r_start, lat_start, lon_start, r, lat, lon );
553  }
554  }
555  }
556  }
557 }
558 
559 
560 
562 
614  Vector& r_v,
615  Vector& lat_v,
616  Vector& lon_v,
617  Vector& za_v,
618  Vector& aa_v,
619  Numeric& lstep,
620  Index& endface,
621  const Numeric& r_start,
622  const Numeric& lat_start,
623  const Numeric& lon_start,
624  const Numeric& za_start,
625  const Numeric& aa_start,
626  const Numeric& ppc,
627  const Numeric& lmax,
628  const Numeric& lat1,
629  const Numeric& lat3,
630  const Numeric& lon5,
631  const Numeric& lon6,
632  const Numeric& r15a,
633  const Numeric& r35a,
634  const Numeric& r36a,
635  const Numeric& r16a,
636  const Numeric& r15b,
637  const Numeric& r35b,
638  const Numeric& r36b,
639  const Numeric& r16b,
640  const Numeric& rsurface15,
641  const Numeric& rsurface35,
642  const Numeric& rsurface36,
643  const Numeric& rsurface16 )
644 {
645  // Radius end latitude of end point
646  Numeric r, lat, lon, l= L_NOT_FOUND; // l not always calculated/needed
647 
648  endface = 0;
649 
650  // Sensor pos and LOS in cartesian coordinates
651  Numeric x, y, z, dx, dy, dz;
652  poslos2cart( x, y, z, dx, dy, dz, r_start, lat_start, lon_start,
653  za_start, aa_start );
654 
655  // Check if crossing with lower pressure level
656  plevel_crossing_3d( r, lat, lon, l, r_start, lat_start, lon_start,
657  za_start, aa_start, x, y, z, dx, dy, dz, ppc,
658  lat1, lat3, lon5, lon6, r15a, r35a, r36a, r16a, true );
659  if( r > 0 )
660  { endface = 2; }
661 
662  // Check if crossing with surface
663  if( rsurface15 >= r15a || rsurface35 >= r35a ||
664  rsurface36 >= r36a || rsurface16 >= r16a )
665  {
666  Numeric rt, latt, lont, lt;
667  plevel_crossing_3d( rt, latt, lont, lt, r_start, lat_start, lon_start,
668  za_start, aa_start, x, y, z, dx, dy, dz, ppc,
669  lat1, lat3, lon5, lon6, rsurface15, rsurface35,
670  rsurface36, rsurface16, true );
671 
672  if( rt > 0 && lt <= l ) // lt<=l to resolve the closest crossing
673  { endface = 7; r = rt; lat = latt; lon = lont; l = lt; }
674  }
675 
676 
677  // Upper pressure level
678  {
679  Numeric rt, latt, lont, lt;
680  plevel_crossing_3d( rt, latt, lont, lt, r_start, lat_start, lon_start,
681  za_start, aa_start, x, y, z, dx, dy, dz, ppc,
682  lat1, lat3, lon5, lon6, r15b, r35b, r36b, r16b,
683  false );
684  if( rt > 0 && lt < l ) // lt<l to resolve the closest crossing
685  { endface = 4; r = rt; lat = latt; lon = lont; l = lt; }
686  }
687 
688  // Latitude 1
689  {
690  Numeric rt, lont, lt;
691  lat_crossing_3d( rt, lont, lt, lat1, lat_start, za_start,
692  x, y, z, dx, dy, dz );
693  if( rt > 0 && lt < l ) // lt<l to resolve the closest crossing
694  { endface = 1; r = rt; lat = lat1; lon = lont; l = lt; }
695  }
696 
697  // Latitude 3
698  {
699  Numeric rt, lont, lt;
700  lat_crossing_3d( rt, lont, lt, lat3, lat_start, za_start,
701  x, y, z, dx, dy, dz );
702  if( rt > 0 && lt < l ) // lt<l to resolve the closest crossing
703  { endface = 3; r = rt; lat = lat3; lon = lont; l = lt; }
704  }
705 
706  // Longitude 5 (only done if solution is lacking)
707  if( !( r>0 && lat>=lat1 && lat<=lat3 && lon>=lon5 && lon<=lon6 ) )
708  {
709  Numeric rt, latt, lt;
710  lon_crossing_3d( rt, latt, lt, lon5, lon_start, za_start, aa_start,
711  x, y, z, dx, dy, dz );
712  if( rt > 0 && lt < l ) // lt<l to resolve the closest crossing
713  { endface = 5; r = rt; lat = latt; lon = lon5; l = lt; }
714  }
715 
716  // Longitude 6 (only done if solution is lacking)
717  if( !( r>0 && lat>=lat1 && lat<=lat3 && lon>=lon5 && lon<=lon6 ) )
718  {
719  Numeric rt, latt, lt;
720  lon_crossing_3d( rt, latt, lt, lon6, lon_start, za_start, aa_start,
721  x, y, z, dx, dy, dz );
722  if( rt > 0 && lt < l ) // lt<l to resolve the closest crossing
723  { endface = 6; r = rt; lat = latt; lon = lon6; l = lt; }
724  }
725 
726  assert( endface );
727 
728  // Check if there is a tangent point inside the grid cell.
729  if( za_start > 90 )
730  {
731  Numeric ltan = geompath_l_at_r( ppc, r_start );
732  if( l-ltan > LACC )
733  {
734  endface = 8;
735  geompath_tanpos_3d( r, lat, lon, l, r_start, lat_start, lon_start,
736  za_start, aa_start, ppc );
737  }
738  }
739 
740  resolve_lon( lon, lon5, lon6 );
741 
742 
743  //--- Create return vectors
744  //
745  Index n = 1;
746  //
747  if( lmax > 0 )
748  {
749  n = Index( ceil( abs( l / lmax ) ) );
750  if( n < 1 )
751  { n = 1; }
752  }
753  //
754  r_v.resize( n+1 );
755  lat_v.resize( n+1 );
756  lon_v.resize( n+1 );
757  za_v.resize( n+1 );
758  aa_v.resize( n+1 );
759  //
760  r_v[0] = r_start;
761  lat_v[0] = lat_start;
762  lon_v[0] = lon_start;
763  za_v[0] = za_start;
764  aa_v[0] = aa_start;
765  //
766  lstep = l / (Numeric)n;
767  //
768  for( Index j=1; j<=n; j++ )
769  {
770  const Numeric lj = lstep * (Numeric)j;
771 
772  cart2poslos( r_v[j], lat_v[j], lon_v[j], za_v[j], aa_v[j],
773  x+dx*lj, y+dy*lj, z+dz*lj, dx, dy, dz, ppc,
774  lat_start, lon_start, za_start, aa_start );
775  resolve_lon( lon_v[j], lon5, lon6 );
776  }
777 
778  //--- Set last point especially, which should improve the accuracy
779  r_v[n] = r;
780  lat_v[n] = lat;
781  lon_v[n] = lon;
782  if( endface == 8 )
783  { za_v[n] = 90; }
784  else
785  { za_v[n] = geompath_za_at_r( ppc, za_start, r_v[n] ); }
786 }
787 
788 
789 
790 
791 
792 
794 
843  Vector& r_v,
844  Vector& lat_v,
845  Vector& za_v,
846  Numeric& lstep,
847  Index& endface,
848  const Numeric& r_start,
849  const Numeric& lat_start,
850  const Numeric& za_start,
851  const Numeric& ppc,
852  const Numeric& lmax,
853  const Numeric& lat1,
854  const Numeric& lat3,
855  const Numeric& r1a,
856  const Numeric& r3a,
857  const Numeric& r3b,
858  const Numeric& r1b,
859  const Numeric& rsurface1,
860  const Numeric& rsurface3 )
861 {
862  // Radius and latitude of end point, and the length to it
863  Numeric r, lat, l= L_NOT_FOUND; // l not always calculated/needed
864 
865  endface = 0;
866 
867  // Check if crossing with lower pressure level
868  plevel_crossing_2d( r, lat, l, r_start, lat_start, za_start, ppc, lat1, lat3,
869  r1a, r3a, true );
870  if( r > 0 )
871  { endface = 2; }
872 
873  // Check if crossing with surface
874  if( rsurface1 >= r1a || rsurface3 >= r3a )
875  {
876  Numeric rt, latt, lt;
877  plevel_crossing_2d( rt, latt, lt, r_start, lat_start, za_start, ppc,
878  lat1, lat3, rsurface1, rsurface3, true );
879 
880  if( rt > 0 && lt <= l ) // lt<=l to resolve the closest crossing
881  { endface = 7; r = rt; lat = latt; l = lt; }
882  }
883 
884  // Upper pressure level
885  {
886  Numeric rt, latt, lt;
887  plevel_crossing_2d( rt, latt, lt, r_start, lat_start, za_start, ppc,
888  lat1, lat3, r1b, r3b, false );
889  if( rt > 0 && lt < l ) // lt<l to resolve the closest crossing
890  { endface = 4; r = rt; lat = latt; /* l = lt; */ }
891  }
892 
893  // Latitude endfaces
894  if( r <= 0 )
895  {
896  if( za_start < 0 )
897  { endface = 1; lat = lat1; }
898  else
899  { endface = 3; lat = lat3; }
900  r = geompath_r_at_lat( ppc, lat_start, za_start, lat );
901  }
902 
903  assert( endface );
904 
905  //Check if there is a tangent point inside the grid cell.
906  bool tanpoint;
907  const Numeric absza = abs( za_start );
908  if( absza > 90 && ( absza - abs(lat_start-lat) ) < 90 )
909  { tanpoint = true; }
910  else
911  { tanpoint = false; }
912 
913  geompath_from_r1_to_r2( r_v, lat_v, za_v, lstep, ppc, r_start, lat_start,
914  za_start, r, tanpoint, lmax );
915 
916  // Force exact values for end point when they are known
917  if( endface == 1 || endface == 3 )
918  { lat_v[lat_v.nelem()-1] = lat; }
919 }
920 
921 
922 
924 
972  Workspace& ws,
973  Array<Numeric>& r_array,
974  Array<Numeric>& lat_array,
975  Array<Numeric>& za_array,
976  Array<Numeric>& l_array,
977  Array<Numeric>& n_array,
978  Array<Numeric>& ng_array,
979  Index& endface,
980  ConstVectorView refellipsoid,
981  ConstVectorView p_grid,
982  ConstVectorView z_field,
983  ConstTensor3View t_field,
984  ConstTensor4View vmr_field,
985  ConstVectorView f_grid,
986  const Numeric& lmax,
987  const Agenda& refr_index_air_agenda,
988  const Numeric& lraytrace,
989  const Numeric& ppc,
990  const Numeric& r_surface,
991  const Numeric& r1,
992  const Numeric& r3,
993  Numeric& r,
994  Numeric& lat,
995  Numeric& za )
996 {
997 // /*
998  // We currently do not handle bending of rays back into the medium due to
999  // refraction. Following check test, whether this is expected to happen for
1000  // the given grid cell and ray zenith angle.
1001  // Generally constant path constant rule applies:
1002  // n1 * r1 * sin(th1) == n3 * r3 * sin(th3),
1003  // where largest path constant values that can occur are (n1*r1) and (n3*r3)
1004  // No solution exists if
1005  // n1 * r1 * sin(th1) > n3 * r3, hence
1006  // n3 < n1 * sin(th1) * r1/r3 (similar for n1 < n3 * sin(th3) * r3/r1)
1007  // Assuming that za is given at r1, we check for n3 >= n1 * sin(za) * r1/r3
1008  Numeric refr1, refr3, refrg;
1009  get_refr_index_1d( ws, refr1, refrg,
1010  refr_index_air_agenda, p_grid, refellipsoid, z_field,
1011  t_field, vmr_field, f_grid, r1 );
1012  get_refr_index_1d( ws, refr3, refrg,
1013  refr_index_air_agenda, p_grid, refellipsoid, z_field,
1014  t_field, vmr_field, f_grid, r3 );
1015  if( refr3 < refr1 * (r1/r3) * sin( DEG2RAD * abs(za) ) )
1016  {
1017  ostringstream os;
1018  os << "For path between r1=" << r1 << "(n-1=" << refr1-1. << ") and r2="
1019  << r3 << "(n-1=" << refr3-1. << "),\n"
1020  << "path calculation will run into refractive back-bending issues.\n"
1021  << "We are currently NOT ABLE to handle them. Consider repeating\n"
1022  << "your calculation without taking refraction into account.";
1023  throw runtime_error(os.str());
1024  }
1025 // */
1026 
1027  // Loop boolean
1028  bool ready = false;
1029 
1030  // Store first point
1031  Numeric refr_index_air, refr_index_air_group;
1032  get_refr_index_1d( ws, refr_index_air, refr_index_air_group,
1033  refr_index_air_agenda, p_grid, refellipsoid, z_field,
1034  t_field, vmr_field, f_grid, r );
1035  r_array.push_back( r );
1036  lat_array.push_back( lat );
1037  za_array.push_back( za );
1038  n_array.push_back( refr_index_air );
1039  ng_array.push_back( refr_index_air_group );
1040 
1041  // Variables for output from do_gridrange_1d
1042  Vector r_v, lat_v, za_v;
1043  Numeric lstep, lcum = 0;
1044 
1045  while( !ready )
1046  {
1047  // Constant for the geometrical step to make
1048  const Numeric ppc_step = geometrical_ppc( r, za );
1049 
1050  // Explicitly check here, that predicted path point is still within
1051  // gridcell and did not undetectedly slip out. This can happen in case of
1052  // close-to-lateral angles and high refraction gradient, which causes a
1053  // bending of the ray back into the medium (kind of reflection). There is
1054  // no proper handling of this back-bending case yet.
1055  if( ( r < r1 - RTOL ) || ( r > r3 + RTOL ) )
1056  {
1057  throw runtime_error(
1058  "Ooops. Path undetectedly left the grid cell.\n"
1059  "This should not happen. But there are issues with cases of high\n"
1060  "refractivity gradients. Seems you unfortunately encountered such\n"
1061  "a case. Little to be done about this now (if this is an option\n"
1062  "for you, run the case without considering refraction). For further\n"
1063  "details, contact Patrick Eriksson.");
1064  }
1065 
1066  // Where will a geometric path exit the grid cell?
1067  do_gridrange_1d( r_v, lat_v, za_v, lstep, endface, r, lat, za, ppc_step,
1068  -1, r1, r3, r_surface );
1069  assert( r_v.nelem() == 2 );
1070 
1071  Numeric za_flagside = za;
1072 
1073  if( lstep <= lraytrace )
1074  {
1075  r = r_v[1];
1076  lat = lat_v[1];
1077  lcum += lstep;
1078  za_flagside = za_v[1];
1079  ready = true;
1080  }
1081  else
1082  {
1083  Numeric l;
1084  if( za <= 90 )
1085  { l = geompath_l_at_r(ppc_step,r) + lraytrace; }
1086  else
1087  {
1088  l = geompath_l_at_r(ppc_step,r) - lraytrace;
1089  if( l < 0 )
1090  { za_flagside = 80; } // Tangent point passed!
1091  }
1092 
1093  r = geompath_r_at_l( ppc_step, l ); // Works als0 for l<0
1094 
1095  lat = geompath_lat_at_za( za, lat,
1096  geompath_za_at_r( ppc_step, za_flagside, r ) );
1097  lcum += lraytrace;
1098  }
1099 
1100  // Refractive index at *r*
1101  get_refr_index_1d( ws, refr_index_air, refr_index_air_group,
1102  refr_index_air_agenda, p_grid, refellipsoid,
1103  z_field, t_field, vmr_field, f_grid, r );
1104 
1105  // Calculate LOS zenith angle at found point.
1106 
1107  const Numeric ppc_local = ppc / refr_index_air;
1108 
1109  if( r >= ppc_local )
1110  { za = geompath_za_at_r( ppc_local, za_flagside, r ); }
1111  else // If moved below tangent point:
1112  {
1113  r = ppc_local;
1114  za = 90;
1115  }
1116 
1117  // Store found point?
1118  if( ready || ( lmax > 0 && lcum + lraytrace > lmax ) )
1119  {
1120  r_array.push_back( r );
1121  lat_array.push_back( lat );
1122  za_array.push_back( za );
1123  n_array.push_back( refr_index_air );
1124  ng_array.push_back( refr_index_air_group );
1125  l_array.push_back( lcum );
1126  lcum = 0;
1127  }
1128  }
1129 }
1130 
1131 
1132 
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
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 za_geom2other_point(const Numeric &r1, const Numeric &lat1, const Numeric &r2, const Numeric &lat2)
za_geom2other_point
The Agenda class.
Definition: agenda_class.h:44
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(Vector &r_v, Vector &lat_v, Vector &za_v, Numeric &lstep, Index &endface, const Numeric &r_start, const Numeric &lat_start, const Numeric &za_start, 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
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
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
void plevel_crossing_3d(Numeric &r, Numeric &lat, Numeric &lon, Numeric &l, const Numeric &r_start0, const Numeric &lat_start, const Numeric &lon_start, const Numeric &za_start, const Numeric &aa_start, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz, const Numeric &ppc, 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 bool &above)
plevel_crossing_3d
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
A constant view of a Tensor4.
Definition: matpackIV.h:141
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:473
const Numeric L_NOT_FOUND
Definition: ppath.cc:93
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
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 cart2poslos(double &r, double &lat, double &lon, double &za, double &aa, const double &x, const double &y, const double &z, const double &dx, const double &dy, const double &dz)
cart2poslos
const Numeric ANGTOL
Definition: ppath.h:103
#define max(a, b)
Definition: continua.cc:20461
void cart2sph(double &r, double &lat, double &lon, const double &x, const double &y, const double &z)
cart2sph
#define abs(x)
Definition: continua.cc:20458
const Numeric R_NOT_FOUND
Definition: ppath.cc:92
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
#define dx
Definition: continua.cc:21928
const Numeric RTOL
Definition: ppath.cc:78
void distance3D(Numeric &l, const Numeric &r1, const Numeric &lat1, const Numeric &lon1, const Numeric &r2, const Numeric &lat2, const Numeric &lon2)
distance3D
Definition: geodetic.cc:658
const Numeric DEG2RAD
Global constant, conversion from degrees to radians.
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 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 refellipsoid, ConstVectorView p_grid, ConstVectorView 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 &ppc, const Numeric &r_surface, const Numeric &r1, const Numeric &r3, Numeric &r, Numeric &lat, Numeric &za)
raytrace_1d_linear_basic
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 lat_crossing_3d(Numeric &r, Numeric &lon, Numeric &l, const Numeric &lat_hit, const Numeric &lat_start, const Numeric &za_start, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz)
lat_crossing_3d
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 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
Workspace class.
Definition: workspace_ng.h:47
void do_gridcell_3d(Vector &r_v, Vector &lat_v, Vector &lon_v, Vector &za_v, Vector &aa_v, Numeric &lstep, Index &endface, const Numeric &r_start, const Numeric &lat_start, const Numeric &lon_start, const Numeric &za_start, const Numeric &aa_start, 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
void lon_crossing_3d(Numeric &r, Numeric &lat, Numeric &l, const Numeric &lon_hit, const Numeric &lon_start, const Numeric &za_start, const Numeric &aa_start, const Numeric &x, const Numeric &y, const Numeric &z, const Numeric &dx, const Numeric &dy, const Numeric &dz)
lon_crossing_3d
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
const Numeric LACC
Definition: ppath.cc:88
void geompath_tanpos_3d(Numeric &r_tan, Numeric &lat_tan, Numeric &lon_tan, Numeric &l_tan, const Numeric &r, const Numeric &lat, const Numeric &lon, const Numeric &za, const Numeric &aa, const Numeric &ppc)
geompath_tanpos_3d
Definition: geodetic.cc:699
void poslos2cart(Numeric &x, Numeric &z, Numeric &dx, Numeric &dz, const Numeric &r, const Numeric &lat, const Numeric &za)
poslos2cart
Definition: geodetic.cc:387
Numeric geompath_l_at_r(const Numeric &ppc, const Numeric &r)
geompath_l_at_r
Definition: ppath.cc:246
const Numeric RAD2DEG
Global constant, conversion from radians to degrees.