利用BASH脚本运行IDL程序、创建子目录或删除文件

有时候你想在一个目录的所有子目录下批量建立文件夹,并且批量的运行这些子目录内新建立好的文件夹内的 IDL 程序,以下是一种可能的解决途径,如果你有更好的方法,欢迎告诉我

Here are the things I want to do:

  1.  I have a data folder containing many sources, like Arp220, NGC981, etc.. Each of them possesses a folder under the whole data folder (./source/Arp220, ./source/NGC891, …). I want to create a subfolder named “plot” under each sources folder, e.g., ./source/Arp220/plot.
  2. And then I want to copy the IDL routine to each “plot” folder I just created and run all the PRO routines at one time.

Here is a bash script being able to do such things:

#!/bin/bash
for i in $(ls -d */); do mkdir -p $i/'plot'; done
for i in $(ls -d */); do cp ./r_name.pro $i/'plot'/; done
for i in $(ls -d */); do cd $PWD/$i/'plot';idl -e ".r ./r_name.pro"; cd ../.. ; done

The first line means creating a “plot” folder under all the subfolders under you present folder. And the next line is to copy the IDL routine (r_name.pro) to each “plot” folder you created. The last line runs all the IDL routines.

顺便再提一个操作,即“批量删除某个目录下所有子目录中包含”X”关键词的文件”。这个可以用一句话可以解决:

find ./* -type f|grep X|while read file; do rm "$file"; done

比如你要删除“./you suck/arXiv.pdf“这种文件。上面那句话中 -type f|grep X|while read file 是为了保证有空格的目录也可以被写入到一个变量里面

在一些大神的提示下,我学到了更简单的做法,比如:

find ./* -type f -name "*X*" -exec rm {} \;

Plotting radial integrated profile from FITS images

Sometimes, you need to know the total power in a particular radius of a FITS image, and the growth of this total power along the radius you choose for checking whether the object is extended or not, or contains some ring or disk structures. Therefore, this script could be very helpful for those kinds of quick checks.

Under the help of Zhangzheng Wen, who gave the idea of calculating the distance from a certain pixel to the image centre, I wrote this IDL radial_power.pro: (Download HERE)

;+
; NAME:
;   RADIAL_POWER
;
; PURPOSE:
;   Calculating the radial_power, returning radius and normalized
;   integrated power within the radius. 
;
; CALLING SEQUENCE:
;   PRO radial_power, image, radius, integrated
;
; INPUTS:
;   image:      read image from a FITS file and put it into a variable.
;
; KEYWORD INPUTS:
;   ...
;
; KEYWORDS SWITCHES:
;   ...
; OUTPUTS:
;   radius:     the radius from center to each circle.
;   integrated: the normalized integrated power within the radius.
;
; EXAMPLES:
;    IDL> image =...
;   ;input the image into the variable after read the image from FITS file or other file format.
;    IDL> radial_power, image, radius, integrated 
;   ;this will return the radius, and integrated.
;
; SIDEEFFECT:
;   ...
;
; RESTRICTIONS:
;   This pro is only useful for one peak structure, namely if
;   you want to study two objects with comparable peak flux, 
;   this pro is not available, but you can edit it to make it
;   function well.
;   NEED "wherenan.pro"
;
; MODIFICATION HISTORY:
;	Write, Chentao YANG (yangcht@gmail.com), 2013-02-04 
;
;-
 
PRO radial_power, image, radius, integrated
 
; make NaN values = 0
IF (SIZE(WHERE(FINITE(image,/NAN)), /N_ELEMENTS) GT 1) THEN image[WHERE(FINITE(image,/NAN))]=0
 
; get the index of the maximum point
s           = SIZE(image,/dimensions)
m           = (WHERE(image EQ MAX(image)))[0]
y_max       = m/s[0]
x_max       = m MOD s[0] 
 
image_size  = SIZE(image)
image_xdist = FINDGEN(image_size[1])#TRANSPOSE(FINDGEN(image_size[2])*0+1)
image_ydist = TRANSPOSE(FINDGEN(image_size[2])#TRANSPOSE(FINDGEN(image_size[1])*0+1))
image_xdist = (image_xdist-x_max)^2
image_ydist = (image_ydist-y_max)^2
image_dist  = SQRT(image_xdist+image_ydist)
 
radius      = FINDGEN(SQRT(MIN([MAX(image_xdist), MAX(image_ydist)]))*0.7)
integrated  = radius
FOR i=0, (SIZE(radius, /N_ELEMENTS)-1) DO BEGIN
  integrated[i]  = TOTAL(image[WHERE(image_dist LE radius[i])])
ENDFOR
integrated  = integrated/MAX(integrated)
 
END

Filling spectral plot (psym=10) in IDL

Two days ago I encountered a problem with plotting spectra using

plot, x, y, PSYM=10

in IDL when I tried to fill the area between the signal and continuum with some color. I firstly tried the POLYFILL, but failed with the PSYM, namely unable to set PSYM=10 if using POLYFILL directly. Then I tried to search if there exists any .pro file written by others for this purpose, but turned out to be nothing. In the end, under the help of Google and Kexin, I found a way through some simple calculation that can solve this by POLYFILL.
Here I attached the .pro file. Hope it can help you as well.

 
;+
; NAME:
;   FILLSPEC
;
; PURPOSE:
;   Two days ago I encountered a problem with plotting spectra using:
;       plot, x, y, PSYM=10
;   in IDL when I tried to fill the area between the signal and continuum 
;   with some color. I firstly tried the POLYFILL, but failed with the PSYM,
;   namely unable to set PSYM=10 if using POLYFILL directly. Then I tried to
;   search if there exists any .pro file written by others for this purpose, 
;   but turned out to be nothing. In the end, under the help of Google and
;   Kexin, I found a way through some simple calculation that can solve this 
;   by POLYFILL.
;
; CALLING SEQUENCE:
;   fillspec, x, y, continuum, psym=, linecolor=, fillcolor=
;
; INPUTS:
;   x, y, continuum. The velocity/frequency/wavelength, flux value and 
;   the continuum value. 
;
; KEYWORD INPUTS:
;   psym      =: Set the plotpsym. For plotting the spectrum, we use 10.
;   linecolor =: Set the spetrum stroke color. 
;   fillcolor =: Set the color that fill the space between flux value
;                and the continuum value.
;
; KEYWORDS SWITCHES:
;   ...
;
; OUTPUTS:
;   ...
;
; EXAMPLES:
;    IDL> x = ...          ; input the velocity/frequency/wavelength array
;    IDL> y = ...          ; input the corresponding flux value array
;    IDL> continuum = ...  ; input the corresponding continuum value (array)
;    IDL> fillspec, x, y, continuum, psym=10, linecolor=220, fillcolor=50
;
; SIDEEFFECT:
;   ...
;
; RESTRICTIONS:
;   (Please check if all the variables are in the same order)
;
; MODIFICATION HISTORY:
;   Write, 2012-04-21, Chentao Yang (yangcht@gmail.com)
;   Only for filling the spectrum after plot the spectrum using 
;   plot and set PSYM as 10.
;   
;   Edited on 2015-01-20, Chentao Yang, 
;   Fix the bug when values beyong YRANGEs are presented. 
;   
;   Edited on 2015-07-17, Chentao Yang, 
;   Fix the bug when the continuum is not a constant value.
;   Fix the bug that this pro may change the variables.
;- 
 
 
PRO fillspec, x0, y0, cont0, LINECOLOR=LINECOLOR, psym=psym, FILLCOLOR=FILLCOLOR
loadct=5
 
x_sort=sort(x0)
x=x0(x_sort)
y=y0(x_sort)
n=n_elements(x)
 
OPLOT, x, y, COLOR=LINECOLOR, psym=psym
 
; Calculate the additional x-values
xnew = (x(0:n-2) + x(1:n-1)) / 2.0D
 
; Create array containing all x-values
xx = REBIN(xnew, N_ELEMENTS(xnew)*2, /sample)
xx = [min(x), min(x), xx, max(x),max(x)]
 
; Create array with the corresponding y-values
yy = REBIN(y, N_ELEMENTS(y)*2, /Sample)
 
; Check if the continuum is a constant
IF (N_ELEMENTS(cont0) NE 1) THEN BEGIN
	cont=cont0(x_sort)
	yy = [cont[0], yy, cont[n_elements(cont)-1]]
ENDIF ELSE BEGIN
	yy = [cont0, yy, cont0]
ENDELSE
 
y_mi_indx = where(yy LT !Y.CRANGE[0])
y_ma_indx = where(yy GT !Y.CRANGE[1])
IF (y_mi_indx[0] GT 0) THEN yy[y_mi_indx]=FLTARR(N_ELEMENTS(y_mi_indx))+!Y.CRANGE[0] 
IF (y_ma_indx[0] GT 0) THEN yy[y_ma_indx]=FLTARR(N_ELEMENTS(y_ma_indx))+!Y.CRANGE[1] 
 
 
 
; Fill the histogram plot
POLYFILL, xx, yy, COLOR=FILLCOLOR
 
; Histogram plot
OPLOT, x, y, COLOR=LINECOLOR, psym=psym
 
END

Here you can download the PRO file.

OR

Clone it.

Least-Squares Fitting of data containing errors

(This is a note personal only! )

There is a powerful tool dealing with data containing both “X” and “Y” errors to make a lease-square fitting in IDL, that is mpfit. This routine allows you to use a self-defined function to do the fitting. Also, you can set up some conditions such as the parameters’ limit or fixed parameters etc. This is fabulous! So beginning with the actual astronomical data have errors both in “X” and “Y”, you need to use a function called “linfitex” (you can see more info from the mpfit website) to both define the function you want to fit and calculate the least-square needed in the fitting. Here is an example of the function:

FUCNTION myfunct, p, x=x, y=y \$
     , sigma_x=sigma_x, sigma_y=sigma_y, _EXTRA=extra
  a = p[0]
  b = p[1]
  f = a*alog10(X)+b
;#########calculate the derivative of the function#########
;##put into resid as (y-f)/sqrt(sigma_y^2+(y'*sigma_x)^2)## 
  resid = (y - f)/sqrt(sigma_y^2 + (a/(x*alog(10))*sigma_x)^2) 
RETURN, resid

In the code above I just defined a function as

    \[y=a \log_{10} (b)+c\]

Then you can read your data. and when you want to fit the data, just using mpfit as:

p      = [1,1,1]
p[1]   = 1         
p[2]   = 2
p[3]   = 3 ;initial guess of the parameters
param  = MPFIT('MYFUNCT', params,\$
      FUNCTARGS={X:x,Y:y,SIGMA_X:x_er,SIGMA_Y:y_er}, PERROR=p_er)

The FUNCTARGS is the place where your data go into the function. And then you can find the results of the parameters from the fitting by print “param” and also the 1σ error from printing p_er. If you only have the error on “Y”, you can use a function like this:

FUNCTION myfunct, p, x=x, y=y, er=er
  a = p[0]     
  b = p[1]
  c = p[2]
  f = a*alog10(b)+c
RETURN, (y-f)/er
END

And then using the same way described above to call “mpfit” and get the results and errors.

Here is also a PYTHON version of this procedure.