187 function etinvert(ich,radi)
200 if (ich < 3 .or. ich > 5)
then
202 .
'-E- : ETINVERT: ICH must be [3,5]. Is ', ich
213 radg = etintegrate(ich,temp);
214 if (radg > radi)
then
218 do while (
abs(radi-radg)/
abs(radi+radg) >= 1e-5)
219 gues = etintegrate(ich,temp+step);
220 if (radg >= radi .and. gues <= radi)
then
221 tnew = temp+step*
abs(radi-radg)/
abs(radg-gues)
223 radg = etintegrate(ich,temp);
226 else if (radg <= radi .and. gues >= radi)
then
227 tnew = temp+step*
abs(radi-radg)/
abs(radg-gues)
229 radg = etintegrate(ich,temp);
237 if (temp < 100.)
then
242 else if (temp > 375.)
then
248 if (radg < radi)
then
262 function etintegrate(ich,etemp)
267 integer,
parameter :: MAX_RSP = 250
277 real*4 resp(max_rsp,3:5);
278 common /resp/ npts, nu0, delnu, resp;
284 real*8 c1,c2,nu,lam,t,w
290 data c1 /1.1910659d-5/
297 w(
tt,nnu) = (c1*nnu*nnu*nnu)/((exp(c2*nnu/
tt)-1.d0))
299 v(
tt,nnu)=sngl((c1*nnu*nnu*nnu*exp(-c2*nnu/
tt))/
300 . (1.d0-(exp(-c2*nnu/
tt))))
302 if (ich < 3 .or. ich > 5)
then
304 .
'-E- : ETINTEGRATE: ICH must be [3,5]. Is ', ich
308 if (npts(ich) > max_rsp)
then
309 write(*,*)
'-E- : ETINTEGRATE: NPTS > MAX_RSP'
312 else if (npts(ich) <= 0)
then
313 write(*,*)
'-E- : ETINTEGEMIS: NPTS <= 0'
318 if (temp < 100.)
then
322 else if (temp > 375.)
then
327 if (nu0(ich) > 100.)
then
329 if (c2*nu/dble(temp) > 20.0d0)
then
332 y(jj) = v(dble(temp),nu) * resp(jj,ich);
333 nu = nu + delnu(ich);
339 y(jj) = sngl(w(dble(temp),nu)) * resp(jj,ich);
340 nu = nu + delnu(ich);
346 lam = dble(nu0(ich));
348 if (c2*nu/dble(temp) > 20.0d0)
then
352 y(jj) = v(dble(temp),nu) * resp(jj,ich);
353 lam = lam + delnu(ich);
360 y(jj) = sngl(w(dble(temp),nu)) * resp(jj,ich);
361 lam = lam + delnu(ich);
370 energy =
simpsn(sngl(delnu(ich)),y,npts(ich),ierr);
371 if (ierr .ne. 0)
then
372 write(*,*)
'-E- : SIMPSN: Integration error..'
376 etintegrate=(energy);
381 subroutine etloadresp(lin, cal)
385 integer,
parameter :: MAX_RSP = 250
393 real*4 resp(max_rsp,3:5);
394 common /resp/ npts, nu0, delnu, resp;
397 character filepath*255;
398 integer*4 ii, jj, ios, ier;
404 character filedir*255
406 data loaded /.false./;
415 call getenv(
'OCDATAROOT', filedir)
416 flen = lenstr(filedir)
417 if (flen .eq. 0)
then
419 .
'-E- ETLOADRESP:Environment variable OCDATAROOT undefined'
424 filepath = filedir(1:flen)//
'/avhrr/cal/'//cal(1:slen)//
'.cal'
425 write(*,*)
'Loading '//filepath(1:lenstr(filepath))
426 open(unit=lin,file=filepath, status=
'OLD');
429 call etgetrsp(lin,ii,npts(ii),resp(1,ii),nu0(ii),
437 factor =
simpsn(sngl(delnu(ii)),resp(1,ii),npts(ii),ier);
439 write(*,*)
'-E- : SIMPSN: Normalizing factor.'
447 resp(jj,ii) = resp(jj,ii) / factor;
458 subroutine etgetrsp(lin,ichn,npts,resp,nu0,delnu,ios)
463 integer,
parameter :: MAX_RSP = 250
482 read(lin,
'(a)',
end=1044) name
483 ln = ftrim(name, 128);
487 else if (name(1:3) ==
"CHN")
then
488 read(name(1:5), 44, iostat=iostat) ic
490 if (ic .eq. ichn)
then
504 read(name(1:ln), 45, iostat=iostat) nu0, delnu, npts
505 if (npts > max_rsp)
then
507 .
'-E- : ETGETRSP: Too many points required (',npts,
508 .
'), maximum is ', max_rsp
515 read(lin,47) (resp(i),i=1,npts)
522 write(*,*)
'-E- : ETGETRSP: No such channel = ',ichn
526 45
format(11x,f14.0,4x,f10.0,5x,i3)
527 47
format(5(e11.5,4x))
532 subroutine etgetvis(lin,cal,slope,intcp,ios)
537 real*4 slope(2), intcp(2);
540 character*(128) name;
541 character filepath*255;
542 character filedir*255;
552 call getenv(
'OCDATAROOT', filedir)
553 flen = lenstr(filedir)
554 if (flen .eq. 0)
then
556 .
'-E- : ETGETVIS: Environment variable OCDATAROOT undefined'
561 filepath = filedir(1:flen)//
'/avhrr/cal/'//cal(1:slen)//
'.cal'
562 write(*,*)
'Loading '//filepath(1:lenstr(filepath))
563 open(unit=lin,file=filepath, status=
'OLD');
565 read(lin,
'(a)',
end=1044) name
566 ln = ftrim(name, 128);
569 else if (name(1:8) ==
"VISIBLE")
then
585 read(lin,43) ic,slope(1),intcp(1);
588 .
'-E- : ETGETVIS: Invalid visible entry for channel 1'
593 read(lin,43) ic,slope(2),intcp(2);
597 .
'-E- : ETGETVIS: Invalid visible entry for channel 2'
609 write(*,*)
'"Visible calibration data not found!'