Blame view

src/idl/dustem_add_linear_params2grid.pro 13 KB
efbec04c   Jean-Philippe Bernard   improved
1
PRO dustem_add_linear_params2grid,input_table_name,pd,iv_min,iv_max,iv_Nvalues,out_filename=out_filename,plog=plog,show_seds=show_seds,help=help
dd80a199   Jean-Philippe Bernard   First commit
2

4225a53d   Jean-Philippe Bernard   improved
3
4
5
6
7
8
9
10
11
12
;+
; NAME:
;    dustem_add_linear_params2grid
; PURPOSE:
;    Adds entries for linear parameters to a Dustemwrap grid fits file
; CATEGORY:
;    Dustem
; CALLING SEQUENCE:
;    dustem_add_linear_params2grid,input_table_name,pd,iv_min,iv_max,iv_Nvalues[,filename=][,plog=][,/show_seds][,/help]
; INPUTS:
907a3ade   Jean-Philippe Bernard   improved
13
;    input_table_name       : input grid fits file name
581fa3e7   Jean-Philippe Bernard   improved
14
15
16
17
;    pd                     : dustemwrap parameter descriptions array
;    iv_min                 : minimum values for parameters in pd
;    iv_max                 : maximum values for parameters in pd
;    iv_Nvalues             : number of parameter values for parameters in pd
4225a53d   Jean-Philippe Bernard   improved
18
; OPTIONAL INPUT PARAMETERS:
581fa3e7   Jean-Philippe Bernard   improved
19
20
;    plog                   : array indicating if the parameters are to be sampled in linear (0) or log scale (1)
;    out_filename           : output fits file name for the table (default ='/tmp/dustem_add_linear_params2grid_final_file.fits')
4225a53d   Jean-Philippe Bernard   improved
21
22
23
24
25
26
27
28
29
30
31
32
33
34
; OUTPUTS:
;    None
; OPTIONAL OUTPUT PARAMETERS:
;    None
; ACCEPTED KEY-WORDS:
;    help      = If set, print this help
; COMMON BLOCKS:
;    None
; SIDE EFFECTS:
;    A file is written
; RESTRICTIONS:
;    The DustEM fortran code must be installed
;    The DustEMWrap IDL code must be installed
; PROCEDURE:
581fa3e7   Jean-Philippe Bernard   improved
35
36
;    Note: A fits file is created at each parameter addition
;    Note: grids do not retain scattering and absorption spectra separately, only totalm extinction. So these are not handled here.
4225a53d   Jean-Philippe Bernard   improved
37
38
39
40
41
42
43
44
45
46
47
48
49
50
; EXAMPLES
;	dustem_init,model='DL07'
;	input_table_name='/Volumes/PILOT_FLIGHT1/PHANGS//ISRF/GRIDS/TEST_DL07_MuseISRF_JWST_G0_4Phangs_noionis_isrfclass15.fits'
;	pd=['(*!dustem_params).grains[0].MDUST_O_MH']
;	default_value=(dustem_get_param_values_default(pd))[0]
;	iv_min=[default_value] & iv_max=[default_value*3.] & iv_Nvalues=[3] & plog=[0]
;	output_table_name='/Volumes/PILOT_FLIGHT1/PHANGS//ISRF/GRIDS/TEST_DL07_MuseISRF_JWST_G0_YPAH_4Phangs_noionis_isrfclass15_essai.fits'
;	dustem_add_linear_params2grid,input_table_name,pd,iv_min,iv_max,iv_Nvalues,filename=output_table_name,plog=plog,/show_seds
; MODIFICATION HISTORY:
;    Written by J.-Ph. Bernard (May 2024)
;    Evolution details on the DustEMWrap gitlab.
;    See http://dustemwrap.irap.omp.eu/ for FAQ and help.  
;-

14188148   Jean-Philippe Bernard   improved
51
52
;JPB note: If fact=1 which happens, we could do nothing (to save time) and just copy that entry of the original table to the new one.
;I started coding this but gave up, as the old seds come out in a format different than needed.
6e36c82b   Jean-Philippe Bernard   improved
53
;JPB note: We could also do without saving intermediate files.
14188148   Jean-Philippe Bernard   improved
54

4225a53d   Jean-Philippe Bernard   improved
55
56
57
58
IF keyword_set(help) THEN BEGIN
  doc_library,'dustem_add_linear_params2grid'
  goto,the_end
ENDIF
dd80a199   Jean-Philippe Bernard   First commit
59

581fa3e7   Jean-Philippe Bernard   improved
60
61
62
63
64
Npd=n_elements(pd)

use_plog=lonarr(Npd)
IF keyword_set(plog) THEN use_plog=plog
;This is the combination values of parameters to be added
dd80a199   Jean-Philippe Bernard   First commit
65
66
parameter_values=dustem_param_range2param_values(iv_min,iv_max,iv_Nvalues,Nc=Nc,log=use_plog)

581fa3e7   Jean-Philippe Bernard   improved
67
68
69
temporary_table_name='/tmp/dustem_add_linear_params2grid_temporary_file.fits'
use_out_filename='/tmp/dustem_add_linear_params2grid_final_file.fits'
IF keyword_set(out_filename) THEN use_out_filename=out_filename
dd80a199   Jean-Philippe Bernard   First commit
70

581fa3e7   Jean-Philippe Bernard   improved
71
72
73
74
IF keyword_set(show_seds) THEN BEGIN
	xrange=[1.,1.e3]
	loadct,13
ENDIF
dd80a199   Jean-Philippe Bernard   First commit
75

581fa3e7   Jean-Philippe Bernard   improved
76
fact_MJy=dustem_get_emission_conversion_factor()
f34148ff   Jean-Philippe Bernard   improved
77
eachi=10L   ;will print progression at eachi
dd80a199   Jean-Philippe Bernard   First commit
78

581fa3e7   Jean-Philippe Bernard   improved
79
80
81
FOR ipd=0L,Npd-1 DO BEGIN
	;==== These are the parameter values of the parameter to be added
	this_parameter_values=dustem_param_range2param_values(iv_min[ipd],iv_max[ipd],iv_Nvalues[ipd],Nc=Nc,log=use_plog[ipd])
e983bf0b   Jean-Philippe Bernard   tried to optimize...
82
83
  use_input_table_name=temporary_table_name
  use_output_table_name=temporary_table_name 	
581fa3e7   Jean-Philippe Bernard   improved
84
85
  IF ipd EQ 0 THEN BEGIN
  	use_input_table_name=input_table_name
e983bf0b   Jean-Philippe Bernard   tried to optimize...
86
87
  ENDIF ;ELSE BEGIN
  IF ipd EQ Npd-1 THEN BEGIN
581fa3e7   Jean-Philippe Bernard   improved
88
	  	use_output_table_name=use_out_filename
581fa3e7   Jean-Philippe Bernard   improved
89
	ENDIF
e983bf0b   Jean-Philippe Bernard   tried to optimize...
90
91
92
93
	;===== read the requested grid fits file
	dustem_read_grid_table2arrays,use_input_table_name $
								  ,Ngrains,model,old_Ncomb,NFparams,INCSP,use_polar,use_double $
								  ,old_parameter_description,old_parameter_pmin,old_parameter_pmax,old_parameter_pnval,old_parameter_plog,old_fparameter_description,old_fparameter_value $
14188148   Jean-Philippe Bernard   improved
94
									,old_param_values,old_Nparams,old_seds,filters,Nfilters,wavs,Nwavs $
cd2aaf3e   Jean-Philippe Bernard   improved
95
								  ,old_I_tot_vec,old_ext_tot_vec,old_I_grains_vec,old_ext_grains_vec,grain_keywords=grain_keywords
e983bf0b   Jean-Philippe Bernard   tried to optimize...
96
	
581fa3e7   Jean-Philippe Bernard   improved
97
98
99
100
101
	new_parameter_description=[old_parameter_description,pd[ipd]]
	new_parameter_pmin=[old_parameter_pmin,iv_min[ipd]]
	new_parameter_pmax=[old_parameter_pmax,iv_max[ipd]]
	new_parameter_pnval=[old_parameter_pnval,iv_Nvalues[ipd]]
	new_parameter_plog=[old_parameter_plog,use_plog[ipd]]
dd80a199   Jean-Philippe Bernard   First commit
102

581fa3e7   Jean-Philippe Bernard   improved
103
104
105
106
  ;This is the new parameter values
  ;It is wrong and used only to save.
  ;This should not be used in that order to write the fits table (see below)
	;new_parameter_values=dustem_param_range2param_values(new_parameter_pmin,new_parameter_pmax,new_parameter_pnval,Nc=Nc,log=new_parameter_plog)
14188148   Jean-Philippe Bernard   improved
107
108
109
  filter_wavs=dustem_filter2wav(filters)
  sed_index=where(filters NE 'SPECTRUM' and filter_wavs gt 0.,count_sed)
  ;stop
dd80a199   Jean-Philippe Bernard   First commit
110

e983bf0b   Jean-Philippe Bernard   tried to optimize...
111
112
	empty_sed=dustem_initialize_sed(Nfilters)
	empty_sed.filter=filters
14188148   Jean-Philippe Bernard   improved
113
	empty_sed.wave=filter_wavs
1f51266a   Jean-Philippe Bernard   removed unnecessa...
114

e983bf0b   Jean-Philippe Bernard   tried to optimize...
115
116
117
118
119
120
  ;==== remove temporary fits file
	fst=file_info(temporary_table_name)
	IF fst.exists EQ 1 THEN BEGIN
		message,'removing '+temporary_table_name,/continue
    str='rm '+temporary_table_name
    spawn,str
9a0f230a   Jean-Philippe Bernard   improved
121
122
123
		message,'removing '+use_out_filename,/continue
    str='rm '+use_out_filename
    spawn,str
e983bf0b   Jean-Philippe Bernard   tried to optimize...
124
125
		;stop
	ENDIF
1f51266a   Jean-Philippe Bernard   removed unnecessa...
126

581fa3e7   Jean-Philippe Bernard   improved
127
128
	New_Ncomb=old_Ncomb*iv_Nvalues[ipd]     ;This is the new number of parameter values combinations
	New_Nparams=old_Nparams+1
dd80a199   Jean-Philippe Bernard   First commit
129

581fa3e7   Jean-Philippe Bernard   improved
130
  ;=== create new arrays for seds, total and partial emission and extinction
934cfd91   Jean-Philippe Bernard   improved custom g...
131
	new_seds_vec=fltarr(Nwavs,New_Ncomb)  ;apparently not used
e983bf0b   Jean-Philippe Bernard   tried to optimize...
132
133
	new_I_tot_vec=fltarr(Nwavs,New_Ncomb)
	new_ext_tot_vec=fltarr(Nwavs,New_Ncomb)
581fa3e7   Jean-Philippe Bernard   improved
134

e983bf0b   Jean-Philippe Bernard   tried to optimize...
135
136
	new_I_grains_vec=fltarr(Nwavs,Ngrains,new_Ncomb)
	new_ext_grains_vec=fltarr(Nwavs,Ngrains,new_Ncomb)
581fa3e7   Jean-Philippe Bernard   improved
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158

	;=== define new quantities
	;=== do the new parameter value structure
	format_template='0.E0'
	IF use_double THEN BEGIN
		format_template='0.D0'
	ENDIF
	str='one_param_st={'
	FOR i=0L,new_Nparams-1 DO BEGIN
	  str=str+'P'+strtrim(i+1,2)+':'+format_template
	  IF i NE new_Nparams-1 THEN str=str+','
	ENDFOR
	str=str+'}'
	toto=execute(str)
	new_param_values=replicate(one_param_st,New_Ncomb)
	;=== fill in the new parameter values structure
	ii=0L
  message,'========= Old parameter values =========',/continue
	FOR i=0L,old_Ncomb-1 DO BEGIN
		print,old_param_values[i]
		FOR j=0L,iv_Nvalues[ipd]-1 DO BEGIN
			FOR k=0L,old_Nparams-1 DO BEGIN
581fa3e7   Jean-Philippe Bernard   improved
159
160
161
162
				new_param_values[ii].(k)=old_param_values[i].(k)   ;This fills in with the content of the fits table
			ENDFOR
			new_param_values[ii].(new_Nparams-1)=(*this_parameter_values[j])[0] ;this adds the new values for the new parameter
			ii=ii+1
dd80a199   Jean-Philippe Bernard   First commit
163
		ENDFOR
dd80a199   Jean-Philippe Bernard   First commit
164
	ENDFOR
581fa3e7   Jean-Philippe Bernard   improved
165
166
  message,'========= New parameter values =========',/continue
  FOR kki=0L,New_Ncomb-1 DO print,new_param_values[kki]
581fa3e7   Jean-Philippe Bernard   improved
167
168
169
170

	new_seds=ptrarr(New_Ncomb)
	em_spectra=ptrarr(New_Ncomb)
	ext_spectra=ptrarr(New_Ncomb)
dd80a199   Jean-Philippe Bernard   First commit
171

581fa3e7   Jean-Philippe Bernard   improved
172
	input_spec=fltarr(Nwavs)   ;input spectrum in MJy/sr
dd80a199   Jean-Philippe Bernard   First commit
173

581fa3e7   Jean-Philippe Bernard   improved
174
175
176
	this_pd=pd[ipd]
	;===== searching for parameter descriptions like '(*!dustem_params).grains[0].MDUST_O_MH'
	pp=strpos(this_pd,'MDUST_O_MH')
dd80a199   Jean-Philippe Bernard   First commit
177

581fa3e7   Jean-Philippe Bernard   improved
178
179
180
181
182
	IF pp NE -1 THEN BEGIN
		;determine which grain type is adressed
		p1=strpos(this_pd,'[')
	  p2=strpos(this_pd,']')
	  igrain=long(strmid(this_pd,p1+1,p2-p1-1))
dd80a199   Jean-Philippe Bernard   First commit
183
184
185
186
		IF incsp NE 1 THEN BEGIN
			message,'Cannot compute new table from grid with no spectra included',/continue
			stop
		ENDIF
581fa3e7   Jean-Philippe Bernard   improved
187
		ii=0L   ;global counter
14188148   Jean-Philippe Bernard   improved
188
		st=dustem_make_empty_fortran_st(Ngrains)
581fa3e7   Jean-Philippe Bernard   improved
189
		;get the default parameter value for the model in use
dd80a199   Jean-Philippe Bernard   First commit
190
191
192
		dustem_init,model=model
		;======= compute the scaling factor
		;== This is the default model value for the parameter to be changed (ie MDUST_O_MH)
581fa3e7   Jean-Philippe Bernard   improved
193
		default_value=(dustem_get_param_values_default(this_pd,param_current_values=param_current_values))[0]
dd80a199   Jean-Philippe Bernard   First commit
194
		FOR i=0L,old_Ncomb-1 DO BEGIN
e983bf0b   Jean-Philippe Bernard   tried to optimize...
195
196
			IF i mod eachi EQ 0 OR i EQ old_Ncomb-1 THEN BEGIN
				message,'Doing initial line '+strtrim(i,2)+' of '+strtrim(old_Ncomb-1,2)+' ('+strtrim(1.*i/(old_Ncomb-1)*100,2)+' %)',/continue
1f51266a   Jean-Philippe Bernard   removed unnecessa...
197
			ENDIF
581fa3e7   Jean-Philippe Bernard   improved
198
			FOR j=0L,iv_Nvalues[ipd]-1 DO BEGIN
e983bf0b   Jean-Philippe Bernard   tried to optimize...
199
				;old_spec=reform(old_I_tot_vec[i,*])
1f51266a   Jean-Philippe Bernard   removed unnecessa...
200
				;stop
581fa3e7   Jean-Philippe Bernard   improved
201
				fact=(*this_parameter_values[j])[0]/default_value    ;Caution, this assumes that the default value was not fixed to a different value when running the grid !
14188148   Jean-Philippe Bernard   improved
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
        ;message,'fact='+strtrim(fact,2),/continue
;stop
        ;IF abs(fact-1.) LE 1.e-5 THEN BEGIN
        ;	message,'Doing nothing ...',/continue
        ;	stop
        ;	new_seds[ii]=ptr_new(old_seds[i])
        ;	FOR kk=0L,Ngrains-1 DO BEGIN
				;		st.sed[*].(kk+1)=old_I_grains_vec[*,kk,ii]/fact_MJy
				;		;st.ext[*].ext_tot=old_ext_grains_vec[*,kk,ii]
				;	ENDFOR
				;	st.sed[*].em_tot=old_I_tot_vec[*,ii]
				;	st.ext[*].ext_tot=old_ext_tot_vec[*,ii]
        ;	em_spectra[ii]=ptr_new(st.sed)
        ;	ext_spectra[ii]=ptr_new(st.ext)
        ;	;stop
        ;	goto,nothing_to_do
        ;ENDIF
e983bf0b   Jean-Philippe Bernard   tried to optimize...
219
220
221
222
223
224
225
226
227
228
229
        ;transfer per-grain emission to new array
				new_I_grains_vec[*,*,ii]=old_I_grains_vec[*,*,i]
				;modifiy per-grain emission for grain type to be modified
				new_I_grains_vec[*,igrain,ii]=old_I_grains_vec[*,igrain,i]*fact
        ;modify total emission spectra accordingly
				new_I_tot_vec[*,ii]=old_I_tot_vec[*,i]-old_I_grains_vec[*,igrain,i]+new_I_grains_vec[*,igrain,ii]
        ;modifiy per-grain extinction for grain type to be modified
				new_ext_grains_vec[*,igrain,ii]=old_ext_grains_vec[*,igrain,i]*fact
				;new_ext_tot_vec[ii,*]=new_ext_tot_vec[ii,*]-old_ext_grains_vec[igrain,i,*]+new_ext_grains_vec[igrain,ii,*]
				new_ext_tot_vec[*,ii]=new_ext_tot_vec[*,ii]-old_ext_grains_vec[*,igrain,i]+new_ext_grains_vec[*,igrain,ii]
				input_spec=new_I_tot_vec[*,ii]
581fa3e7   Jean-Philippe Bernard   improved
230
				;=== Note: The st structure is useful to construct, even if not used in dustem_compute_sed, to get the right shape for the spectra
14188148   Jean-Philippe Bernard   improved
231
232
233
				;IF ii EQ 0 THEN BEGIN
				;	st=dustem_make_empty_fortran_st(Ngrains)
				;ENDIF
581fa3e7   Jean-Philippe Bernard   improved
234
				;=== fill in st with the new emission and extinction spectrum
e983bf0b   Jean-Philippe Bernard   tried to optimize...
235
				;=== reform is needed here because of the structure aspect of st
581fa3e7   Jean-Philippe Bernard   improved
236
				FOR kk=0L,Ngrains-1 DO BEGIN
e983bf0b   Jean-Philippe Bernard   tried to optimize...
237
					st.sed[*].(kk+1)=new_I_grains_vec[*,kk,ii]/fact_MJy
14188148   Jean-Philippe Bernard   improved
238
					;st.ext[*].ext_tot=new_ext_grains_vec[*,kk,ii]
dd80a199   Jean-Philippe Bernard   First commit
239
				ENDFOR
1f51266a   Jean-Philippe Bernard   removed unnecessa...
240
241
242
				st.sed[*].em_tot=st.sed[*].(1) ;This is the emission of grain type 1
				FOR kk=1L,Ngrains-1 DO BEGIN
					st.sed[*].em_tot=st.sed[*].em_tot+st.sed[*].(kk+1)
dd80a199   Jean-Philippe Bernard   First commit
243
				ENDFOR
14188148   Jean-Philippe Bernard   improved
244
				st.ext[*].ext_tot=new_ext_tot_vec[*,ii]
581fa3e7   Jean-Philippe Bernard   improved
245
				;=== This below is why it is useful to get st. Caution, st.sed must be in fortran units, not MJy/sr
dd80a199   Jean-Philippe Bernard   First commit
246
247
				em_spectra[ii]=ptr_new(st.sed)
				ext_spectra[ii]=ptr_new(st.ext)
581fa3e7   Jean-Philippe Bernard   improved
248
				;At that stage, st.sed.em_tot should equal input_spec, which is not the case
e983bf0b   Jean-Philippe Bernard   tried to optimize...
249
250
				;diff_spec=input_spec-st.sed.em_tot*fact_MJy
				;dustem_Ised=dustem_compute_sed([0.],st=st,filters=filters)   ;finally recompute the SED (should give the same result as below)
581fa3e7   Jean-Philippe Bernard   improved
251
				;The line below gives the same result as the above, which is a confirmation that st and input_spec are consistent
14188148   Jean-Philippe Bernard   improved
252
				dustem_Ised=dustem_compute_sed([0.],input_spec=input_spec,filters=filters,sed_index=sed_index,/no_sort)   ;finally recompute the SED
e983bf0b   Jean-Philippe Bernard   tried to optimize...
253
				new_sed=empty_sed
dd80a199   Jean-Philippe Bernard   First commit
254
				new_sed.stokesI=dustem_Ised
e983bf0b   Jean-Philippe Bernard   tried to optimize...
255
				new_sed.sigmaII=0.01*dustem_Ised
dd80a199   Jean-Philippe Bernard   First commit
256
				new_seds[ii]=ptr_new(new_sed)
14188148   Jean-Philippe Bernard   improved
257
258
259

				nothing_to_do:

dd80a199   Jean-Philippe Bernard   First commit
260
				IF keyword_set(show_seds) THEN BEGIN
581fa3e7   Jean-Philippe Bernard   improved
261
				  colors=long(range_gen(New_Ncomb,[0,255]))
dd80a199   Jean-Philippe Bernard   First commit
262
263
264
265
266
267
					wavs=(*new_seds[ii]).wave
					Is=(*new_seds[ii]).stokesI
	    			IF ii EQ 0 THEN BEGIN
					    tit='Grid SEDs model '+model
					    xtit='Wavelength [mic]'
					    ytit='SED [MJy/sr for NH=1.e20 H/cm^2]'
dd80a199   Jean-Philippe Bernard   First commit
268
269
270
271
272
				    	cgplot,wavs,Is,psym=-4,/ylog,yr=yrange,xtit=xtit,ytit=ytit,tit=tit,/xlog,xrange=xrange,/xsty
			    	ENDIF ELSE BEGIN
			    		cgoplot,wavs,Is,psym=-4,color=colors[ii]
			    	ENDELSE
			    	ishow=100
dd80a199   Jean-Philippe Bernard   First commit
273
274
275
276
				ENDIF
				ii=ii+1
			ENDFOR
		ENDFOR
581fa3e7   Jean-Philippe Bernard   improved
277
278
279
280
281
282
	ENDIF ELSE BEGIN
		message,'parameter description '+this_pd,' not processed',/continue
		stop
	ENDELSE
	;===== write the output grid
	;help,model,new_parameter_description,new_parameter_pmin,new_parameter_pmax,new_parameter_pnval,new_parameter_plog,new_param_values,new_seds
dd80a199   Jean-Philippe Bernard   First commit
283

581fa3e7   Jean-Philippe Bernard   improved
284
285
  ;recompute the pointer to new_parameter_values from effective parameter values
  ;This is needed to have the grid saved actually be saved with the right parameter order ...
f34148ff   Jean-Philippe Bernard   improved
286
  ;message,'========= New parameter values written to table =========',/continue
581fa3e7   Jean-Philippe Bernard   improved
287
288
289
290
291
292
293
294
  new_parameter_values=ptrarr(New_Ncomb)
  new_Nparams=old_Nparams+1
  vec=dblarr(new_Nparams)
  FOR ii=0L,New_Ncomb-1 DO BEGIN
  	FOR kk=0L,new_Nparams-1 DO BEGIN
  		vec[kk]=new_param_values[ii].(kk)
  	ENDFOR
  	new_parameter_values[ii]=ptr_new(vec)
1f51266a   Jean-Philippe Bernard   removed unnecessa...
295
  	;print,*new_parameter_values[ii]
581fa3e7   Jean-Philippe Bernard   improved
296
297
  ENDFOR
 	;stop
dd80a199   Jean-Philippe Bernard   First commit
298

e983bf0b   Jean-Philippe Bernard   tried to optimize...
299
300
  message,'Writing '+use_output_table_name,/continue

581fa3e7   Jean-Philippe Bernard   improved
301
302
303
304
305
306
307
308
309
	dustem_write_sed_grid,model,new_parameter_description, $
						 new_parameter_pmin,new_parameter_pmax,new_parameter_pnval,new_parameter_plog, $
						 new_parameter_values,new_seds, $
						 em_spectra=em_spectra, $
						 ext_spectra=ext_spectra, $
						 fpd=old_fparameter_description, $
						 fiv=old_fparameter_value, $
						 filename=use_output_table_name, $
						 use_polarisation=use_polarisation, $
cd2aaf3e   Jean-Philippe Bernard   improved
310
311
						 double=use_double, $
						 grain_keywords=grain_keywords
581fa3e7   Jean-Philippe Bernard   improved
312
313
	 ;stop
ENDFOR
dd80a199   Jean-Philippe Bernard   First commit
314

4225a53d   Jean-Philippe Bernard   improved
315
the_end:
dd80a199   Jean-Philippe Bernard   First commit
316
317

END