2011年7月27日水曜日

2次元 SMM fit

#!/Library/Frameworks/EPD64.framework/Versions/Current/bin/python

SCRIPT_TYPE="InstrumentalityOfMankind" 
SCRIPT_SHELL="python" 
SCRIPT_NAME="Fit_infile_.py"
SCRIPT_HELP_ARG=""
SCRIPT_HELP_SENTENCE=""
SCRIPT_NUM_ARG=0
SCRIPT_VERSION=1.0

###IMPORT###
import SMM
from scipy.optimize import leastsq 
from numpy import *
import sys

###HELP###
if len(sys.argv[1:])!=SCRIPT_NUM_ARG:
    print 'Name: '+SCRIPT_NAME
    print 'Arguments: '+SCRIPT_HELP_ARG
    print 'Explanation: '+SCRIPT_HELP_SENTENCE
    sys.exit()


###MAIN###
def SMM_norm_free(norm_nsc,norm_nsd,norm_gb,norm_gd):
    return lambda l,b: norm_nsc * SMM.SMM(l,b,'NSC') + norm_nsd * SMM.SMM(l,b,'NSD') +  norm_gb * SMM.SMM(l,b,'GB') + norm_gd * SMM.SMM(l,b,'GD')

def SMM_norm_free_array(param,l_array,b_array):
    res=[]
    num=0
    for l in l_array :
        b=b_array[num]
        res.append(SMM_norm_free(*param)(l,b))
        num+=1                
    return res

def fit_SMM_norm_free(param,l_array,b_array,sb_array,sb_err_array):
    errorfunction = lambda p: (array(SMM_norm_free_array(p,l_array,b_array))-array(sb_array))/sb_err_array
    p=leastsq(errorfunction,param,full_output=True)
    return p

def demo(para=[1.0,2.0,3.0,4.0],para0=[1.1,1.2,1.3,1.4]) :
    l_in=linspace(0.0,5.0,10)
    b_in=linspace(0.0,5.0,10)
    data = SMM_norm_free_array(para,l_in, b_in)
    data_err=0.001*array(data)
    param=fit_SMM_norm_free(para0,l_in,b_in,data,data_err)
    
    val=param[0]
    covar=param[1]
    val_err=[]
    for num in range(0,4):
        val_err.append((covar[num][num])**0.5)
    
    print val,val_err
    

def fit_SMM_indata(infile):
    para0=(1.0,1.0,1.0,1.0)
    data=loadtxt(infile)
    l_in=data[:,0]
    b_in=data[:,1]
    sb_in=data[:,2]
    sb_in_err=data[:,3]
    param=fit_SMM_norm_free(para0,l_in,b_in,sb_in,sb_in_err)
    
    
    val=param[0]
    covar=param[1]
    val_err=[]
    for num in range(0,4):
        val_err.append((covar[num][num])**0.5)
    
    print val,val_err



fit_SMM_indata(str(sys.argv[1]))