subroutine fetch_from_best_fit_model - (pixel_num,pol_ch_num,frequency, - file_of_model_sets, - use_only_latest_model, - model_val,scaled_model_val, - fractional_error, model_was_used) implicit none real*4 frequency, - model_val, - scaled_model_val, - fractional_error integer*4 pixel_num, - pol_ch_num logical*4 use_only_latest_model, - model_was_used character*(*) file_of_model_sets c+ c_name fetch_from_best_fit_model c_function return the the best-fit-model-based value and its scaled version c (e.g. Tcal (in K) values of both high and low CAL) c for a requested ALFA pixel (0-6), polarization channel (1 or 2) and c the frequency in MHz. c Data consisting of one or several sets of best-fit coefficients c describing multi-order polymonial & harmonic fits to the c relevant measurement results for several bands, c dual linear polarization channels and the seven ALFA pixels c are read from a disk-file (named in "file_of_model_sets"). c These data are used to compute the values of the fitted parameter c (and its scaled version), and when more than one measurements are available, c a suitably weighted average estimate is computed. If the model is unavailable for c the sought parameter combination, a default value is returned. c c_call call fetch_from_best_fit_model(pixel_num,pol_ch_num,frequency, c_call - file_of_model_sets, c_call - use_only_latest_model, c_call - model_val,scaled_model_val, c_call - fractional_error, model_was_used) c_/pixel_num i*4 ALFA pixel number (0-6) c_/pol_ch_num i*4 polarization channels number (1 or 2, for pol A & B respectively) c_/frequency r*4 frequency in (MHz) c_/file_of_model_sets ch*(*) path/name for the disk-file containing sets of best-fit model_coefficients c_/use_only_latest_model l*4 logical variable that is set to .true. to force usage of only the latest model c_/model_val r*4 returned value for primary (modeled) parameter (e.g. Tcal-high) c_/scaled_model_val r*4 returned value for scaled version of the modeled parameter (e.g. Tcal-low) c_/fractional_error r*4 returned estimate of fractional error in the parameters (both, primary & scaled) c_/model_was_used l*4 returned logical flag to signal if the returned value is indeed based on the model c_ i.e. .false. would indicate that the returned value is the default one. c_author desh c_created 11-SEP-2004 c_latest 11-SEP-2004 c_reference Arecibo Observatory, prompted by ALFA c_keyword system calibration, Cal noise temperature, spectral scans, polynomial/harmonic fit c- c local variables integer*4 maxpixels, maxpols, maxbands, maxpar parameter (maxpixels = 8, maxpols = 2, - maxbands = 10, maxpar = 128 ) integer*4 number, ipix, istk, i, iparam, iset, n_summed, - mm, nn, start_set, loglun/51/, - order_p, order_h, in_pix, nchar, - no_of_sets(maxpixels,maxpols), - polynomial_order(maxpixels,maxpols,maxbands), - harmonic_order(maxpixels,maxpols,maxbands) real*4 param_read(maxpixels,maxpols,maxbands,maxpar), - freq_start(maxpixels,maxpols,maxbands), - freq_end(maxpixels,maxpols,maxbands), - ref_freq(maxpixels,maxpols,maxbands), - freq_scale_factor(maxpixels,maxpols,maxbands), - rms_uncertainty(maxpixels,maxpols,maxbands), - weightage(maxbands), weight_sum, cal_sum, freq_arg, - rms_now(maxbands), cal_now(maxbands), - atemp, btemp, rms, f_offset, f_factor, - param(maxpar), scaling_factor, additional_offset, - default_value, pi, pi_by_2, degrees_of_freedom logical*4 data_already_loaded/.false./, skip character*128 already_loaded_file/'NOT LOADED'/ character*20 data_source_tag save no_of_sets, polynomial_order, harmonic_order, - param_read, freq_start, freq_end, ref_freq, - freq_scale_factor, rms_uncertainty, - data_already_loaded, default_value, - already_loaded_file, scaling_factor, - pi_by_2, pi,additional_offset pi = acos(-1.) pi_by_2 = pi/2. in_pix = pixel_num + 1 ! conversion to ensure Fortran-compatible indexing (i.e. starting from 1) if(data_already_loaded)then ! assess if the requested file is different from the one already read nn = nchar(already_loaded_file) mm = nchar(file_of_model_sets) if(nn.ne.mm)then ! the filenames differ at the basic level data_already_loaded = .false. ! we need to open the new file else ! check to see if the names differ, though the names are equal in length if(file_of_model_sets(1:mm).ne. - already_loaded_file(1:nn))then ! the names do differ, if not their lengths data_already_loaded = .false. ! we need to open the new file end if end if end if if(.not.data_already_loaded)then ! try to open the file and read all the fit-parameter sets i = -1 open(unit=loglun, - file=file_of_model_sets(1:nchar(file_of_model_sets)), - type='old',err=100) i = 1 100 if(i.lt.0)then ! failed to open the file; report it and stop write(*,*)'ERROR in opening ', - file_of_model_sets(1:nchar(file_of_model_sets)) stop ! no point in continuing else ! file is open; read it now i = -1 read(loglun,*,err=102,end=102)default_value,additional_offset, - scaling_factor i = 1 102 if(i.lt.0)then ! filed to read the first line write(*,*)'ERROR in reading 1st line from ', - file_of_model_sets(1:nchar(file_of_model_sets)) stop ! no point in continuing end if c first initialize the band-counter array do ipix=1,maxpixels do istk=1,maxpols no_of_sets(ipix,istk) = 0 end do end do do while(.true.) ! read till the expected format is conformed to read(loglun,*,err=101,end=101) - ipix,istk,atemp,btemp, - f_offset,f_factor, - order_p,order_h,mm, - (param(nn),nn=1,mm), - rms, degrees_of_freedom - ,data_source_tag(1:4) ipix = ipix + 1 ! conversion to ensure Fortran-compatible indexing no_of_sets(ipix,istk) = no_of_sets(ipix,istk) + 1 iset = no_of_sets(ipix,istk) freq_start(ipix,istk,iset) = atemp freq_end(ipix,istk,iset) = btemp ref_freq(ipix,istk,iset) = f_offset freq_scale_factor(ipix,istk,iset) = f_factor c note that the total number of coefficients mm = order_p + 2*order_h, c and that order_p accounts for the zeroth-order coefficient polynomial_order(ipix,istk,iset) = order_p - 1 harmonic_order(ipix,istk,iset) = order_h do nn=1,mm param_read(ipix,istk,iset,nn) = param(nn) end do rms_uncertainty(ipix,istk,iset) = rms end do 101 close(unit=loglun) write(33,*)'Successfully read data from ', - file_of_model_sets(1:nchar(file_of_model_sets)) do ipix=1,maxpixels do istk=1,maxpols write(33,*)'read ',ipix,istk,no_of_sets(ipix,istk) cc write(33,*)'read ',ipix,istk,no_of_sets(ipix,istk) end do end do data_already_loaded = .true. ! note the new status ! and the name of the loaded file already_loaded_file(1:) = - file_of_model_sets(1:nchar(file_of_model_sets)) end if end if c Initialise some arrays/variables used in the fitting routines if(no_of_sets(in_pix,pol_ch_num).gt.0)then ! data are readily available model_was_used = .true. ! let us begin with some optimism if(use_only_latest_model)then start_set = no_of_sets(in_pix,pol_ch_num) else start_set = 1 end if n_summed = 0 do iset=start_set,no_of_sets(in_pix,pol_ch_num) skip = .false. if((frequency.ge.freq_start(in_pix,pol_ch_num,iset)).and. - (frequency.le.freq_end(in_pix,pol_ch_num,iset)))then if(model_was_used)n_summed = n_summed + 1 ! otherwise overwrite the earlier one that was not within the model range c perform scaling of the frequency argument for range -1 to + 1 freq_arg = - (frequency - ref_freq(in_pix,pol_ch_num,iset)) - /freq_scale_factor(in_pix,pol_ch_num,iset) rms_now(n_summed) = rms_uncertainty(in_pix,pol_ch_num,iset) c see if the rms is sensible. If not, force a large value to it if(rms_now(n_summed).le.0.0)then rms_now(n_summed) = default_value + additional_offset end if model_was_used = .true. else ! frequency is outside the available range, appeal for estimate at the nearest freq. if(n_summed.eq.0)then ! arrange to pick some value at the nearest-freq., if no good ones are yet found n_summed = 1 if(frequency.le.freq_start(in_pix,pol_ch_num,iset))then freq_arg = (freq_start(in_pix,pol_ch_num,iset) - - ref_freq(in_pix,pol_ch_num,iset)) - /freq_scale_factor(in_pix,pol_ch_num,iset) else freq_arg = (freq_end(in_pix,pol_ch_num,iset) - - ref_freq(in_pix,pol_ch_num,iset)) - /freq_scale_factor(in_pix,pol_ch_num,iset) end if c admit that uncertainity may be large in the blind assignment rms_now(n_summed) = default_value + additional_offset c and deny the optimism model_was_used = .false. else skip = .true. end if end if if(.not.skip)then c nominal weight corresponding to inverse of the associated noise variance weightage(n_summed) = - 1./(rms_now(n_summed)*rms_now(n_summed)) c compute this estimate from the best-fit co-efficients c starting with the zero-th order component atemp = param_read(in_pix,pol_ch_num,iset,1) iparam = 1 c then, add the first and the higher-order polynomial terms if(polynomial_order(in_pix,pol_ch_num,iset).gt.0)then do i = 1,polynomial_order(in_pix,pol_ch_num,iset) iparam = iparam + 1 atemp = atemp - + param_read(in_pix,pol_ch_num,iset,iparam) - * freq_arg**real(i) end do end if c finally, include also the harmonic-pair contribution if(harmonic_order(in_pix,pol_ch_num,iset).gt.0)then do i = 1,harmonic_order(in_pix,pol_ch_num,iset) c first cos term iparam = iparam + 1 atemp = atemp - + param_read(in_pix,pol_ch_num,iset,iparam) - *cos(real(i)*pi_by_2*freq_arg) c and then the sin term iparam = iparam + 1 atemp = atemp - + param_read(in_pix,pol_ch_num,iset,iparam) - *sin(real(i)*pi_by_2*freq_arg) end do end if cal_now(n_summed) = atemp end if end do if(.not.model_was_used)then ! freq seems to be outside the available range, and would have n_summed.eq.1 write(33,*)'No data available in this frequency range !' model_val = cal_now(1) + additional_offset scaled_model_val = model_val*scaling_factor fractional_error = 1.0 ! one may specify another number that may be meaningful, but we will assume 100% error model_was_used = .false. return else ! we do have some data if(n_summed.eq.1)then model_val = cal_now(1) + additional_offset scaled_model_val = model_val*scaling_factor fractional_error = rms_now(1)/model_val else ! more than one measurements available weight_sum = 0.0 cal_sum = 0.0 btemp = 0.0 do i=1,n_summed cal_sum = cal_sum + cal_now(i)*weightage(i) weight_sum = weight_sum + weightage(i) btemp = btemp + (rms_now(i)*weightage(i))**2 end do if(weight_sum.gt.0.0)then model_val = cal_sum/weight_sum + additional_offset fractional_error = sqrt(btemp)/(weight_sum*model_val) else model_val = default_value + additional_offset fractional_error = 1.0 ! one may specify another number that may be meaningful, but we will assume 100% error end if scaled_model_val = model_val*scaling_factor end if end if else ! necessary data are not readily available write(33,*)'No data available for this pixel/pol_channel ', - pixel_num,pol_ch_num model_val = default_value + additional_offset scaled_model_val = default_value*scaling_factor fractional_error = 1.0 ! one may specify another number that may be meaningful, but we will assume 100% error model_was_used = .false. return end if return end c============================================================================================== cccc Integer function NCHAR(string) ccccC ccccC Routine to count the number of characters in the ccccC input string. Looks for the last occurrence of ccccC non-(null, blank or tab character) ccccC ccccC ccccC cccc Implicit none cccc integer*4 i,ipos cccc character*(*) string cccc character blank,tab,null,c ccccC data blank,tab,null/' ',9,0/ cccc blank=' ' cccc tab=char(9) cccc null=char(0) cccc cccc ipos = 0 cccc i = len(string) cccc do while (i.gt.0.and.ipos.eq.0) cccc c = string(i:i) cccc if (c.ne.blank.and.c.ne.tab.and.c.ne.null) ipos = i cccc i = i - 1 cccc end do cccc cccc NCHAR = ipos cccc return cccc end c==============================================================================================