OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
earth.f
Go to the documentation of this file.
1  subroutine earth( pos, vel, widphse1, widphfl1, widphse2,
2  1 widphfl2, yawest, navctl, nad_bod, nadbodfl, ssdec)
3 c
4 c earth( pos, vel, widphse1, widphfl1, widphse2,
5 c widphfl2, yawest, navctl, nad_bod, nadbodfl)
6 c
7 c purpose: determine the nadir body vectors from the earth
8 c sensor data
9 c
10 c calling arguments:
11 c
12 c name Type i/o description
13 c -------- ---- --- -----------
14 c pos r*4 i size 3 orbit position vector(km)
15 c vel r*4 i size 3 orbit velocity vector(km/sec)
16 c widphse1 r*4 i size 2 processed earth width and
17 c phase for sensor 1, in degrees
18 c widphfl1 i*4 i flag for earth sensor 1 - 0 good
19 c widphse2 r*4 i size 2 processed earth width and
20 c phase for sensor 2, in degrees
21 c widphfl2 i*4 i flag for earth sensor 2 - 0 good
22 c yawest r*4 i estimated spacecraft yaw in degrees
23 c navctl struct i controls including specification of
24 c cone angle, crossing biases, matrix
25 c for transform of sensor frame to
26 c spacecraft frame
27 c nad_bod r*4 o size 3 unit nadir vector in spacecraft frame
28 c nadbodfl i*4 o flag for the the nadir vector, 1 - bad
29 c ssdec r*4 o sine of the solar declination
30 c
31 c by: w. robinson, gsc, 12 apr 93 from version by j whiting
32 c
33 c notes: this was taken from the original earth without much mods.
34 c due to the flags for earth sensors and deterministic nadir
35 c solution added, it may be good to re-write this from scratch
36 c wr
37 c
38 c modification history:
39 c w. robinson, gsc, 16 apr 93 this is a change in body of code
40 c from original: flip indicies on array hv so that vector
41 c routines may be called with this array.
42 c f. patt, gsc, 23 mar 95: fixed scanner azimuth calculation to
43 c allow for arbitrary scanner orientation
44 c f. patt, gsc, 23 sep 97: added horizon sensor in- and out-crossing
45 c angle calibrations as 4th-order polynomials
46 c f. patt, gsc, 16 oct 97: corrected a bug in the oblateness modeling.
47 c f. patt, saic gsc, 24 apr 98: added horizon crossing biases to
48 c initial calculation of rho at scanner axis azimuth
49 c f. patt, saic gsc, 11 may 98: modified calculation of horizon angle
50 c by adjusting flattening factor to fit apparent horizon
51 c scanner triggering height.
52 c f. patt, saic gsc, 19 june 98: further modified calculation of the
53 c horizon angle to include a correction based on the solar
54 c declination angle and the subsatellite latitude.
55 c f. patt, saic gsc, 5 feb 99: set nad_bod to zero if flags are set.
56 c f. patt, saic gsc, 23 aug 1999: modified horizon sensor calibrations
57 c to use vendor-provided calibration tables; modified horizon
58 c angle calculation seasonal correction based on analysis of
59 c sensor data with revised calibration.
60 c f. patt, saic gsc, 9 dec 1999: modified horizon angle calculation
61 c seasonal correction based on further analysis of sensor data.
62 c f. patt, saic gsc, 20 dec 1999: modified flattening factor used for
63 c horizon angle calculation based on analysis of island targets.
64 c f. patt, saic gsc, 7 jan 2000: incorporated new form of horizon angle
65 c calculation to allow seasonal shifting of ellipsoid in place
66 c of adjustment of subsatellite radius.
67 c f. patt, saic gac, 4 apr 2001: modified flattening factor to include
68 c seasonal correction based on island target results.
69 
70  implicit none
71 c
72 #include "nav_cnst.fin"
73 #include "navctl_s.fin"
74 c
75  type(navctl_struct) :: navctl
76 c
77  integer*4 nadbodfl, widphfl1, widphfl2
78  real*4 pos(3),vel(3),yawest,nad_bod(3),
79  1 widphse1(2), widphse2(2)
80 c
81  real*8 xnad(3)
82 
83  real*4 xang(4), v2(3), xang1(4)
84  real*4 sca1, sca2, cca1, cca2, hv0(3), hv(3,4)
85  real*4 sc0x(3), sc0y(3), sc0z(3), sc0zm, pxy2, posm
86  real*4 n(3), nm, ndoty, ndotz, yawref
87  real*4 aa1, aa2, at1, da1, da2, at2, rh1, rh2, xaz(4)
88  real*4 lat, clat, ffrom1, ffrom2, clat2, re2
89  real*4 slat, cxaz, cxaz2, rho(4), crho
90  real*4 rho1, rho2, hv1(3), hv2(3), or1(3), or2(3), or3(3),
91  1 norm1, norm2, norm3, vmag, nad1(3), ff, ssdec, ssdfac(3)
92  real*4 aa, ab, ac, af, ah, ai, aj, dx, ftmp
93  real*8 mv(3), state(3,3), xnadm
94  real*8 a, b, c
95  integer*2 i,j,k, fl1, fl2, iax
96  integer*4 ierr,i3
97 c
98  data rho/4*0.0/
99  data ssdfac/0.0,11.5,-4.5/
100  data i3/3/
101 c
102 c
103 c
104 c start, set local flags from the input flags and see
105 c if this index should be done
106 c
107  fl1 = widphfl1
108  fl2 = widphfl2
109 c
110  if( ( fl1 .eq. 1 ) .and. ( fl2 .eq. 1 ) ) then
111  nadbodfl = 1
112  do i=1,3
113  nad_bod(i) = 0.0
114  end do
115  else
116  nadbodfl = 0
117 
118 c initialize:
119 c convert scanner alignment angles to transformation matrices
120 c
121 c convert earth width and split-to-index to in- &
122 c out-crossing angles(xang)
123 c
124  if( fl1 .eq. 0 ) then
125  xang1(1) = widphse1(2) - widphse1(1)/2.d0
126  xang1(2) = widphse1(2) + widphse1(1)/2.d0
127  else
128  xang1(1) = -135.
129  xang1(2) = -135.
130  end if
131 c
132  if( fl2 .eq. 0 ) then
133  xang1(3) = widphse2(2) - widphse2(1)/2.d0
134  xang1(4) = widphse2(2) + widphse2(1)/2.d0
135  else
136  xang1(3) = 135.
137  xang1(4) = 135.
138  end if
139 
140 c apply horizon angle calibrations
141 
142  call earcal(xang1, xang)
143 
144  do i=1,4
145  xang(i) = xang(i)/radeg
146  xang1(i) = xang1(i)/radeg
147  end do
148 
149 c compute horizon vectors in scanner coordinates
150 c and immediately transform to spacecraft coordinates
151 c hv = [cos(xang)sin(ca), sin(xang)sin(ca), cos(ca)]
152 c
153  if( fl1 .eq. 0 ) then
154  sca1 = sin( navctl%ear1sca / radeg )
155  cca1 = cos( navctl%ear1sca / radeg )
156  do i=1,2
157  hv0(1) = cos(xang(i))*sca1
158  hv0(2) = sin(xang(i))*sca1
159  hv0(3) = cca1
160  call matvec( navctl%ear_mat(1,1,1), hv0, hv(1,i) )
161  enddo
162  end if
163 c
164 c note that we define both scanner systems approx paallel
165 c within the limits of the orientation
166 c
167  if( fl2 .eq. 0 ) then
168  sca2 = sin( navctl%ear2sca / radeg )
169  cca2 = cos( navctl%ear2sca / radeg )
170  do i=3,4
171  hv0(1) = cos(xang(i))*sca2
172  hv0(2) = sin(xang(i))*sca2
173  hv0(3) = cca2
174  call matvec( navctl%ear_mat(1,1,2), hv0, hv(1,i) )
175  enddo
176  end if
177 
178 c compute model nadir-to-horizon angles(rho)
179 c corresponding to vectors:
180 c compute(zero attitude) s/c axes from pos. & vel. vectors
181 c
182  v2(1) = vel(1) - omegae*pos(2)
183  v2(2) = vel(2) + omegae*pos(1)
184  v2(3) = vel(3)
185 
186  call crossp(pos,v2,sc0z)
187  sc0zm = sqrt(sc0z(1)*sc0z(1)+sc0z(2)*sc0z(2)+sc0z(3)*sc0z(3))
188  pxy2 = pos(1)*pos(1)+pos(2)*pos(2)
189  posm = sqrt(pxy2+pos(3)*pos(3))
190  do i=1,3
191  sc0x(i) = -pos(i)/posm
192  sc0z(i) = sc0z(i)/sc0zm
193  enddo
194  call crossp(sc0z,sc0x,sc0y)
195 
196 c compute local north vector(=<pos>x<earth centered z-axis>x<pos>)
197 c
198  n(1) = -pos(1)*pos(3)
199  n(2) = -pos(2)*pos(3)
200  n(3) = pxy2
201  nm = sqrt(n(1)*n(1)+n(2)*n(2)+n(3)*n(3))
202  do i=1,3
203  n(i) = n(i)/nm
204  enddo
205 
206 c compute s/c(geocentric) latitude
207 c
208  lat = asin(pos(3)/posm)
209 
210 c compute sub-satellite earth radius(wertz p102, eq 4-23)
211 c
212  clat = cos(lat)
213 
214 c flattening factor with seasonal correction
215  ff = 1.d0/1.85d2 - ssdec/6.4d2
216 c ff = 1.d0/2.0d2
217  ffrom1 = 1.d0 - ff
218  ffrom2 = 2.d0 - ff
219  ftmp = ffrom2*ff/(ffrom1*ffrom1)
220  clat2 = clat*clat
221 
222 c set-up to compute rho values(wertz p102, eq 4-24)
223 c
224  re2 = (re + ssdfac(3)*ssdec)**2
225  slat = pos(3)/posm
226  aa = 1.0
227  ab = 1.0 + ftmp*clat*clat
228  ac = 1.0 + ftmp*slat*slat
229  af = 2.*ftmp*slat*clat
230  ah = posm*af
231  ai = 2.*posm*ac
232  aj = posm*posm*ac - re2
233 
234 c add seasonal correction to ellipsoid
235 
236  dx = ssdfac(1) + ssdfac(2)*ssdec
237  ah = ah - 2.*dx*clat/(ffrom1*ffrom1)
238  ai = ai - 2.*dx*slat/(ffrom1*ffrom1)
239  aj = aj - (2.*dx*pos(3) - dx*dx)/(ffrom1*ffrom1)
240 c compute yaw reference azimuth(angle betw. north & orbit y-axis)
241 c
242  ndoty = n(1)*sc0y(1)+n(2)*sc0y(2)+n(3)*sc0y(3)
243  ndotz = n(1)*sc0z(1)+n(2)*sc0z(2)+n(3)*sc0z(3)
244  yawref = -atan2(ndotz,ndoty)
245 
246 c compute scanner axis azimuth(scanner axis align.:
247 c yaw~= +90 & +270)
248 c
249  da1 = atan2(navctl%ear_mat(3,3,1),navctl%ear_mat(2,3,1))
250  aa1 = yawref + yawest/radeg + da1
251  da2 = atan2(navctl%ear_mat(3,3,2),navctl%ear_mat(2,3,2))
252  aa2 = yawref + yawest/radeg + da2
253 
254 c compute in- & out-xing azimuths(xaz) relative to north
255 c xaz = aa +/- atan(tan(ca)sin(ew/2)) (approx)
256 c
257  if( fl1 .eq. 0 ) then
258 
259 c first compute rho at scanner axis azimuth
260  cxaz = cos(aa1)
261  cxaz2 = cxaz*cxaz
262  a = ai*ai - 4.d0*aj*ac
263  b = (4.d0*aj*af - 2.d0*ah*ai)*cxaz
264  c = (ah*ah+4.d0*aj*(aa-ab))*cxaz2 - 4.d0*aj*aa
265  rh1 = atan(2.d0*a/(sqrt(b*b - 4.d0*a*c) - b))
266 c need to include average crossing bias
267  rh1 = rh1 + (navctl%e1biasic + navctl%e1biasoc)/(2.d0*radeg)
268 
269 c compute azimuth offsets to horizon crossings
270  at1 = acos((sca1**2*cos(xang(2)-xang(1))
271  * + cca1**2-cos(rh1)**2) / sin(rh1)**2)/2.
272  xaz(1) = aa1 + at1
273  xaz(2) = aa1 - at1
274  end if
275 c
276  if( fl2 .eq. 0 ) then
277 
278 c first compute rho at scanner axis azimuth
279  cxaz = cos(aa2)
280  cxaz2 = cxaz*cxaz
281  a = ai*ai - 4.d0*aj*ac
282  b = (4.d0*aj*af - 2.d0*ah*ai)*cxaz
283  c = (ah*ah+4.d0*aj*(aa-ab))*cxaz2 - 4.d0*aj*aa
284  rh2 = atan(2.d0*a/(sqrt(b*b - 4.d0*a*c) - b))
285 c need to include average crossing bias
286  rh2 = rh2 + (navctl%e2biasic + navctl%e2biasoc)/(2.d0*radeg)
287 
288 c compute azimuth offsets to horizon crossings
289  at2 = acos((sca2**2*cos(xang(4)-xang(3))
290  * + cca2**2 - cos(rh2)**2) / sin(rh2)**2)/2.
291  xaz(3) = aa2 + at2
292  xaz(4) = aa2 - at2
293  end if
294 
295 c compute nadir-to-horizon angles(rho):
296 c
297  if( fl1 .eq. 0 ) then
298  do i=1,2
299  cxaz = cos(xaz(i))
300  cxaz2 = cxaz*cxaz
301  a = ai*ai - 4.d0*aj*ac
302  b = (4.d0*aj*af - 2.d0*ah*ai)*cxaz
303  c = (ah*ah+4.d0*aj*(aa-ab))*cxaz2 - 4.d0*aj*aa
304  rho(i) = atan(2.d0*a/(sqrt(b*b - 4.d0*a*c) - b))
305  enddo
306  end if
307 c
308  if( fl2 .eq. 0 ) then
309  do i=3,4
310  cxaz = cos(xaz(i))
311  cxaz2 = cxaz*cxaz
312  a = ai*ai - 4.d0*aj*ac
313  b = (4.d0*aj*af - 2.d0*ah*ai)*cxaz
314  c = (ah*ah+4.d0*aj*(aa-ab))*cxaz2 - 4.d0*aj*aa
315  rho(i) = atan(2.d0*a/(sqrt(b*b - 4.d0*a*c) - b))
316  enddo
317  end if
318 
319 c add horizon scanner triggering biases to angles
320 
321  rho(1) = rho(1) + navctl%e1biasic/radeg
322  rho(2) = rho(2) + navctl%e1biasoc/radeg
323  rho(3) = rho(3) + navctl%e2biasic/radeg
324  rho(4) = rho(4) + navctl%e2biasoc/radeg
325 
326  if (navctl%lvdbug.gt.2) then
327  do i=1,4
328  write (68,*) (hv(j,i),j=1,3),rho(i), xang(i)
329  end do
330  end if
331 
332 c depending on whether one or 2 earth sensors are
333 c available, process differently
334 c
335  if( fl1 .eq. 1 .or. fl2 .eq. 1 ) then
336 c
337 c use deterministic approach. first, move the good
338 c horizon and rho vectors to new locations
339 c
340  if( fl1 .eq. 0 ) then
341 c
342  do iax = 1, 3
343  hv1(iax) = hv(iax,1)
344  hv2(iax) = hv(iax,2)
345  end do
346  rho1 = rho(1)
347  rho2 = rho(2)
348 c
349  else
350 c
351  do iax = 1, 3
352  hv1(iax) = hv(iax,3)
353  hv2(iax) = hv(iax,4)
354  end do
355  rho1 = rho(3)
356  rho2 = rho(4)
357 c
358  end if
359 c
360 c construct an orthonormal system out of the
361 c horizon vectors: sum, difference and cross product
362 c
363  do iax = 1, 3
364  or1(iax) = hv1(iax) + hv2(iax)
365  or3(iax) = hv1(iax) - hv2(iax)
366  end do
367 c
368  call crossp( hv1, hv2, or2 )
369 c
370  norm1 = vmag( or1 )
371  norm2 = vmag( or2 )
372  norm3 = vmag( or3 )
373 c
374  do iax = 1, 3
375  or1(iax) = or1(iax) / norm1
376  or2(iax) = or2(iax) / norm2
377  or3(iax) = or3(iax) / norm3
378  end do
379 c
380 c dot the nadir vector into the or1 and or3 axis
381 c to get the components of the nadir in this system
382 c( remember, cos(rho) = nadir * horiz vec )
383 c
384  nad1(1) = ( cos( rho1 ) + cos( rho2 ) ) / norm1
385  nad1(2) = ( cos( rho1 ) - cos( rho2 ) ) / norm3
386  nad1(3) = -sqrt( 1 - nad1(1) * nad1(1) - nad1(2) * nad1(2) )
387 c
388 c the square root has 2 solutions. pick the one that has
389 c the greatest x component: thus, if or2(1) is < 0,
390 c use the negative solution
391 c
392 c if( or2(1) .lt. 0 ) nad1(3) = -nad1(3)
393 c
394 c transform the nadir vector back to spacecraft system
395 c
396  do iax = 1, 3
397  nad_bod(iax) = nad1(1) * or1(iax) +
398  1 nad1(2) * or3(iax) +
399  1 nad1(3) * or2(iax)
400  end do
401 c
402 c
403  else
404 c
405 c in this case, three or four valid vectors, use
406 c least-squares algorithm:
407 c initialize(upper right of) state matrix &
408 c measurement vector as zeroes
409 c
410  do i=1,3
411  mv(i)=0.d0
412  do j=i,3
413  state(i,j)=0.d0
414  enddo
415  enddo
416 
417 c for each(valid) horizon vector:
418 c
419  do i=1,4
420 
421 c update state matrix(upper right half) using
422 c horizon vector(hv(i)):
423 c [sm] = [sm] + <hv>#<HV> (# denotes the outer product)
424 c
425  do j=1,3
426  do k=j,3
427  state(j,k) = state(j,k) + hv(j,i)*hv(k,i)
428  enddo
429  enddo
430 
431 c update measurement vector using horizon
432 c vector(hv) & angle(rho):
433 c <mv> = <mv> + cos(rho)*<hv>
434 
435  crho = cos(rho(i))
436  do j=1,3
437  mv(j) = mv(j) + crho*hv(j,i)
438  enddo
439  enddo
440 
441 c update lower left hand of(symmetric) state matrix
442 c
443  state(2,1) = state(1,2)
444  state(3,1) = state(1,3)
445  state(3,2) = state(2,3)
446 
447 c invert state matrix and solve state equations for nadir vector:
448 c <nad> = inv[sm]*<mv>
449 c
450  call invert(state,mv,i3,i3,xnad,ierr)
451 
452 c normalize the nadir vector to be safe.
453 c
454  xnadm = dsqrt( xnad(1) * xnad(1) +
455  1 xnad(2) * xnad(2) +
456  1 xnad(3) * xnad(3))
457  do i=1,3
458  nad_bod(i) = xnad(i)/xnadm
459  enddo
460 c
461  end if
462  end if
463 c
464  return
465  end
466 
467  subroutine earcal(xang,yang)
468 
469 c
470 c earcal( xang, yang)
471 c
472 c purpose: apply calibrations to the horizon scanner crossing angles
473 c
474 c calling arguments:
475 c
476 c name Type i/o description
477 c -------- ---- --- -----------
478 c xang(4) r*4 i input horizon scanner crossing angles
479 c(order is hsa in, out, hsb in, out
480 c yang(4) r*4 o output(calibrated) angles
481 c
482 c by: f. s. patt, saic gsc, may 4, 1998
483 c
484 c notes: this subroutine uses the look-up tables provided by osc (and,
485 c presumably, the manufacturer).
486 c
487 c modification history:
488 c
489 c modified to use cubic lagrange instead of linear interpolation.
490 c f. s. patt, saic gsc, august 16, 1999
491 c
492 c modified to read all calibration tables from a single file.
493 c f. s. patt, saic gsc, august 18, 1999
494 
495  implicit none
496 
497  real*4 xang(4), yang(4)
498  real*4 hscal(541,4), hsref(541), xfac, xf(4)
499  integer*4 i, j, lenstr
500  character*256 hsfile
501  logical first
502 
503  common /calcom/ hscal, hsref, first
504 
505  data first/.true./
506 
507 c If first call, open files and read tables
508 
509 
510  if (first) then
511  first = .false.
512  hsfile = '$NAVCTL/'
513  call filenv(hsfile, hsfile)
514  hsfile = hsfile(1:lenstr(hsfile)) // 'hs_cal.dat'
515  open(1,file = hsfile)
516  do j=1,541
517  read(1,*) hsref(j),(hscal(j,i), i=1,4)
518  end do
519 
520  close(1)
521  end if
522 
523 c find location in lookup table for each angle
524  do i=1,4
525  j = 3
526  if (xang(i).lt.0) xang(i) = xang(i) + 360.
527  dowhile((xang(i).gt.hsref(j)).and.(j.lt.540))
528  j = j + 1
529  end do
530 
531 c compute interpolation factors
532  xfac = (xang(i)-hsref(j-1))/(hsref(j)-hsref(j-1))
533  xf(1) = -xfac*(xfac-1.0)*(xfac-2.0)/6.
534  xf(2) = (xfac+1.0)*(xfac-1.0)*(xfac-2.0)/2.
535  xf(3) = -(xfac+1.0)*xfac*(xfac-2.0)/2.
536  xf(4) = (xfac+1.0)*xfac*(xfac-1.0)/6.
537 
538  yang(i) = xf(1)*hscal(j-2,i) + xf(2)*hscal(j-1,i)
539  * + xf(3)*hscal(j,i) + xf(4)*hscal(j+1,i)
540 
541  end do
542 
543  return
544  end
subroutine earcal(xang, yang)
Definition: earth.f:468
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame index(0-based) was previously being used the frame number(1-based) should have been used. There were only a few minor changes to the uncertainty index(maximum of 1 digit). 3. Some inefficient arrays(Sigma_RVS_norm_sq) were eliminated and some code lines in Preprocess_L1A_Data were moved into Process_OBCEng_Emiss. There were no changes to the product. Required RAM was reduced by 20 MB. Now
float * vector(long nl, long nh)
Definition: nrutil.c:15
===========================================================================V4.1.3 12/18/2002============================================================================Changes which do not affect scientific output:1. The R *LUT was eliminated and the equivalent formulation for R *, i.e. 1/(m1 *e_sun_over_pi), was substituted for it in the only instance of its use, which is in the calculation of the RSB uncertainty index. This reduces the size of the Reflective LUT HDF file by approximately 1/4 to 1/3. The equivalent formulation of R *differed from the new by at most 0.056% in test granules and uncertainty differences of at most 1 count(out of a range of 0-15) were found in no more than 1 in 100, 000 pixels. 2. In Preprocess.c, a small error where the trailing dropped scan counter was incremented when the leading dropped scan counter should have been was fixed. This counter is internal only and is not yet used for any purpose. 3. NEW MYD02OBC Metadata Configuration Files. MCST wishes to have the OBC files archived even when the Orbit Number is recorded as "-1". Accordingly, ECS has delivered new MCF files for OBC output having all elements in the OrbitCalculatedSpatialDomain container set to "MANDATORY=FALSE". 4. pgs_in.version is now reset to "1" in Metadata.c before the call to look up the geolocation gringpoint information.============================================================================V4.1.1 CODE SPECIFIC TO MODIS/AQUA(FM1) 10/03/2002============================================================================Two changes were made to the code which do not affect scientific output:1. A bug which caused PGE02 to fail when scans were dropped between granules was fixed.(The length of the error message generated was shortened.) 2. Messages regarding an invalid MCST LUT Version or an invalid Write High Resolution Night Mode Output value in the PCF file were added.==============================================================================V4.1.0 CODE SPECIFIC TO MODIS/AQUA(FM1)(NEVER USED IN PRODUCTION) 07/30/2002==============================================================================Changes which impact scientific output of code:1. The LUT type of the RVS corrections was changed to piecewise linear. In addition the RVS LUTs were changed from listing the RVS corrections to listing the quadratic coefficients necessary to make the RVS corrections. The coefficients are now calculated by interpolating on the granule collection time and the RVS corrections are then generated using the interpolated coefficients. Previously used Emissive and Reflective RVS LUT tables were eliminated and new ones introduced. Several changes were made to the code which should not affect scientific output. They are:1. The ADC correction algorithm and related LUTs were stripped from the code.(The ADC correction has always been set to "0" so this has no scientific impact.) 2. Some small changes to the code, chiefly to casting of variables, were added to make it LINUX-compatible. Output of code run on LINUX machines displays differences of at most 1 scaled integer(SI) from output of code run on IRIX machines. The data type of the LUT "dn_sat_ev" was changed to float64 to avoid discrepancies seen between MOD_PR02 run on LINUX systems and IRIX systems where values were flagged under one operating system but not the other. 3. Checking for non-functioning detectors, sector rotation, incalculable values of the Emissive calibration factor "b1", and incalculable values of SV or BB averages was moved outside the loop over frames in Emissive_Cal.c since none of these quantities are frame-dependent. 4. The code was altered so that if up to five scans are dropped between the leading/middle or middle/trailing granules, the leading or trailing granule will still be used in emissive calibration to form a cross-granule average. QA bits 25 and 26 are set for a gap between the leading/middle and middle/trailing granules respectively. This may in rare instances lead to a change in emissive calibration coefficients for scans at the beginning or end of a granule. 5.(MODIS/AQUA ONLY) The name of the seed(error message) file was changed from "MODIS_36100.h" to "MODIS_36110.h". 6. Metadata.c was changed so that the source of the geolocation metadata is the input geolocation file rather than the L1A granule. 7. To reduce to overall size of the reflective LUT HDF files, fill values were eliminated from all LUTs previously dimensioned "BDSM"([NUM_REFLECTIVE_BANDS] *[MAX_DETECTORS_PER_BAND] *[MAX_SAMPLES_PER_BAND] *[NUM_MIRROR_SIDES]) in the LUT HDF files. Each table piece is stored in the HDF file with dimensions NUM_REFLECTIVE_INDICES, where NUM_REFLECTIVE_INDICES=[NUM_250M_BANDS *DETECTORS_PER_250M_BAND *SAMPLES_PER_250M_BAND *NUM_MIRROR_SIDES]+[NUM_500M_BANDS *DETECTORS_PER_500M_BAND *SAMPLES_PER_500M_BAND *NUM_MIRROR_SIDES]+[NUM_1000M_BANDS *DETECTORS_PER_1KM_BAND *SAMPLES_PER_1KM_BAND *NUM_MIRROR_SIDES] with SAMPLES_PER_250M_BAND=4, SAMPLES_PER_500M_BAND=2, and SAMPLES_PER_1KM_BAND=1. Values within each table piece appear in the order listed above. The overall dimensions of time dependent BDSM LUTs are now[NUM_TIMES] *[NUM_REFLECTIVE_INDICES], where NUM_TIMES is the number of time dependent table pieces. 8. Checking for non-functioning detectors, sector rotation, incalculable values of the Emissive calibration factor "b1", and incalculable values of SV or BB averages was moved outside the loop over frames in Emissive_Cal.c since none of these quantities are frame-dependent. 9. The code was altered so that if up to five scans are dropped between the leading/middle or middle/trailing granules, the leading or trailing granule will still be used in emissive calibration to form a cross-granule average. QA bits 25 and 26 are set for a gap between the leading/middle and middle/trailing granules respectively. This may in rare instances lead to a change in emissive calibration coefficients for scans at the beginning or end of a granule. 10. The array of b1s in Preprocess.c was being initialized to -1 outside the loop over bands, which meant that if b1 could not be computed, the value of b1 from the previous band for that scan/detector combination was used. The initialization was moved inside the band loop.============================================================================V3.1.0(Original Aqua-specific code version) 02/06/2002============================================================================AQUA-Specific changes made:1. A correction to a problem with blackbody warmup on bands 33, 35, and 36 was inserted. PC Bands 33, 35, and 36 on MODIS Aqua saturate on BB warmup before 310K, which means current code will not provide correct b1 calibration coefficients when the BB temperatures are above the saturation threshold. A LUT with default b1s and band-dependent temperature thresholds will be inserted in code. If the BB temperature is over the saturation threshold for the band, the default b1 from the table is used. 2. The number of possible wavelengths in the Emissive LUT RSR file was changed to 67 in order to accommodate the Aqua RSR tables. 3. Several changes to the upper and lower bound limits on LUT values were inserted. Changes to both Aqua and Terra Code:1. A check was put into Emissive_Cal.c to see whether the value of b1 being used to calibrate a pixel is negative. If so, the pixel is flagged with the newly created flag TEB_B1_NOT_CALCULATED, value 65526, and the number of pixels for which this occurs is counted in the QA_common table. 2. The array of b1s in Preprocess.c was being initialized to -1 outside the loop over bands, which meant that if b1 could not be computed, the value of b1 from the previous band for that scan/detector combination was used. The initialization was moved inside the band loop. 3. Minor code changes were made to eliminate compiler warnings when the code is compiled in 64-bit mode. 4. Temperature equations were upgraded to be MODIS/Aqua or MODIS/Terra specific and temperature conversion coefficients for Aqua were inserted.========================================================================================================================================================ALL CHANGES BELOW ARE TO COMMON TERRA/AQUA CODE USED BEFORE 02/06/2002========================================================================================================================================================v3.0.1 11/26/2001============================================================================Several small changes to the code were made, none of which changes the scientific output:1. The code was changed so that production of 250m and 500m resolution data when all scans of a granule are in night mode may be turned off/on through the PCF file. 2. A check on the times of the leading and trailing granules was inserted. If a leading or trailing granule does not immediately precede or follow(respectively) the middle granule, it is treated as a missing granule and a warning message is printed. 3. The code now reads the "MCST Version Number"(e.g. "3.0.1.0_Terra") from the PCF file and checks it against the MCST Version number contained in the LUT HDF files. This was done to allow the user to make sure the code is being run using the correct LUT files.(The designators "0_Terra", "1_Terra", etc.) refer to the LUT versions.) 4. A small bug in Preprocess.c was corrected code
Definition: HISTORY.txt:661
void remember(char *Vdata_name, int32 Vdata_id)
Definition: remember.c:8
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT table(e.g. m1) to be used for interpolation. The table values will be linearly interpolated using the tables corresponding to the node times bracketing the granule time. If the granule time falls before the time of the first node or after the time of the last node
subroutine invert(a, b, n, l, c, ier)
Definition: invert.f:3
#define flattening
Definition: vincenty.c:24
void initialize(int pixref_flag, int blkref_flag)
Definition: Usds.c:1371
subroutine earth(pos, vel, widphse1, widphfl1, widphse2,
Definition: earth.f:2
float ** matrix(long nrl, long nrh, long ncl, long nch)
Definition: nrutil.c:60
void fit(float x[], float y[], int ndata, float sig[], int mwt, float *a, float *b, float *siga, float *sigb, float *chi2, float *q)
#define real
Definition: DbAlgOcean.cpp:26
Definition: novas.h:66
subroutine filenv(infil, outfil)
Definition: filenv.f:2
real *4 function vmag(vec)
Definition: vmag.f:2
int state(double tjdTDB, JPLIntUtilType *util, double posvel[13][6], double *pnut)
subroutine matvec(xm, v1, v2)
Definition: matvec.f:2
#define re
Definition: l1_czcs_hdf.c:701
===========================================================================V5.0.48(Terra) 03/20/2015 Changes shown below are differences from MOD_PR02 V5.0.46(Terra)============================================================================Changes noted for V6.1.20(Terra) below were also instituted for this version.============================================================================V6.1.20(Terra) 03/12/2015 Changes shown below are differences from MOD_PR02 V6.1.18(Terra)============================================================================Changes from v6.1.18 which may affect scientific output:A situation can occur in which a scan which contains sector rotated data has a telemetry value indicating the completeness of the sector rotation. This issue is caused by the timing of the instrument command to perform the sector rotation and the recording of the telemetry point that reports the status of sector rotation. In this case a scan is considered valid by L1B and pass through the calibration - reporting extremely high radiances. Operationally the TEB calibration uses a 40 scan average coefficient, so the 20 scans(one mirror side) after the sector rotation are contaminated with anomalously high radiance values. A similar timing issue appeared before the sector rotation was fixed in V6.1.2. Our analysis indicates the ‘SET_FR_ENC_DELTA’ telemetry correlates well with the sector rotation encoder position. The use of this telemetry point to determine scans that are sector rotated should fix the anomaly occured before and after the sector rotation(usually due to the lunar roll maneuver). The fix related to the sector rotation in V6.1.2 is removed in this version.============================================================================V6.1.18(Terra) 10/01/2014 Changes shown below are differences from MOD_PR02 V6.1.16(Terra)============================================================================Added doi attributes to NRT(Near-Real-Time) product.============================================================================V6.1.16(Terra) 01/27/2014 Changes shown below are differences from MOD_PR02 V6.1.14(Terra)============================================================================Migrate to SDP Toolkit 5.2.17============================================================================V6.1.14(Terra) 06/26/2012 Changes shown below are differences from MOD_PR02 V6.1.12(Terra)============================================================================Added the doi metadata to L1B product============================================================================V6.1.12(Terra) 04/25/2011 Changes shown below are differences from MOD_PR02 V6.1.8(Terra)============================================================================1. The algorithm to calculate uncertainties for reflective solar bands(RSB) is updated. The current uncertainty in L1B code includes 9 terms from prelaunch analysis. The new algorithm regroups them with the new added contributions into 5 terms:u1:the common term(AOI and time independent) and
Definition: HISTORY.txt:126
subroutine crossp(v1, v2, v3)
Definition: crossp.f:2
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DWITH_MPI") target_link_libraries(afrt_nc4 $
Definition: CMakeLists.txt:16
subroutine single
Definition: single.f:2
#define eq(a, b)
Definition: rice.h:167
subroutine phase
Definition: phase.f:2
flags
Definition: DDAlgorithm.h:22
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be calibrated(i.e., the BB DNs are not used as a backup for the SWIR bands). 3. Piecewise linear LUT capability was added to the code. If a LUT table is marked as piecewise linear
subroutine linear(xp, x, y, n, yp)
Definition: afrt_rt1.f:1255
#define f
Definition: l1_czcs_hdf.c:702
void radius(double A)
Definition: proj_report.c:132
#define degrees(radians)
Definition: niwa_iop.c:30
this program makes no use of any feature of the SDP Toolkit that could generate such a then geolocation is calculated at that and then aggregated up to Resolved feature request Bug by adding three new int8 SDSs for each high resolution offsets between the high resolution geolocation and a bi linear interpolation extrapolation of the positions This can be used to reconstruct the high resolution geolocation Resolved Bug by delaying cumulation of gflags until after validation of derived products Resolved Bug by setting Latitude and Longitude to the correct fill resolving to support Near Real Time because they may be unnecessary if use of entrained ephemeris and attitude data is turned on(as it will be in Near-Real-Time processing).
algorithm
Definition: DDProcess.h:25