dustem_first_guess.pro
2.37 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
FUNCTION dustem_first_guess,pd,pval,sed,pd_to_update,pd_filter_names, $
fpd=fpd, $
fiv=fpval, $
pol=pol, $
verbose=verbose, $
use_previous_fortran=use_previous_fortran, $
dustem_predicted_sed=dustem_predicted_sed, $
output_dustem_st=output_dustem_st, $
input_dustem_st=input_dustem_st, $
no_reset_plugin_structure=no_reset_plugin_structure
;computes new parameter initial values, based on provided SED and provided filter names associated to parameters
;where a linear calculation can be done
;sed is actually not used. Uses sed in !dustem_data instead.
Npd=n_elements(pd_to_update)
Nfilters=n_elements(pd_filter_names)
IF Npd NE Nfilters THEN BEGIN
message,'pd_to_update and pd_filter_names must be of same length',/continue
stop
ENDIF
new_pval=pval ;defaults is no change to provideid initial values
;==== compute SED on pd_filter_names with provided pd,iv
;dummy=0 ;otherwise model is not recomputed !!
IF not keyword_set(input_dustem_st) THEN st=0 ELSE st=input_dustem_st
dustem_init_params,!dustem_model,pd,pval,fpd=fpd,fiv=fpval,pol=pol,no_reset_plugin_structure=no_reset_plugin_structure
dustem_Ised=dustem_compute_sed(pval,st=st,use_previous_fortran=use_previous_fortran)
output_dustem_st=st ;This is to output the Fortran structure
;=== for keyword_output. Caution, this is the SED for the input parameter values, not the output ones
dustem_predicted_sed=dustem_Ised
;=== The above two are not used
;sed_filters=(*(*!dustem_data).sed).filt_names
;sed_values=(*(*!dustem_data).sed).filt_names
;Caution: below, dustem_Ised and sed.filter can have different size, leading to an error,
;in particular when use_previous_fortran, due to filtering by dustem_set_data. Below is a (very dirty) fix
data_filter_names=(*(*!dustem_data).sed).filt_names ;CAUTION can be different from sed.filter
use_sed=sed
;use_sed=(*(*!dustem_data).sed).values
FOR i=0L,Nfilters-1 DO BEGIN
;dirty fix
;ind1=where(sed.filter EQ pd_filter_names[i],count1)
ind1=where(data_filter_names EQ pd_filter_names[i],count1)
;stop
IF count1 NE 0 THEN BEGIN
;rap=dustem_Ised[ind1[0]]/sed[ind1[0]].stokesI
rap=dustem_Ised[ind1[0]]/use_sed[ind1[0]].stokesI
ind2=where(pd EQ pd_to_update[i],count2)
new_pval[ind2[0]]=pval[ind2[0]]/rap
ENDIF
ENDFOR
IF keyword_set(verbose) THEN BEGIN
print,pdverbose
print,pval
print,new_pval
ENDIF
;stop
RETURN,new_pval
END