<s> import sys , getopt , os <EOL> sys . path . append ( '<STR_LIT>' ) <EOL> import numpy as np <EOL> import matplotlib as mpl <EOL> import matplotlib . pyplot as plt <EOL> import pyfits <EOL> import exorings <EOL> import j1407 <EOL> from scipy . optimize import fmin <EOL> from scipy . interpolate import UnivariateSpline <EOL> from scipy . interpolate import interp1d <EOL> from matplotlib . patches import PathPatch <EOL> mpl . interactive ( True ) <EOL> mpl . rc ( '<STR_LIT:image>' , interpolation = '<STR_LIT>' , origin = '<STR_LIT>' , cmap = '<STR_LIT>' ) <EOL> mpl . rc ( '<STR_LIT>' , limits = ( - <NUM_LIT:7> , <NUM_LIT:7> ) ) <EOL> G = <NUM_LIT> <EOL> yr = <NUM_LIT> * <NUM_LIT> <EOL> msol = <NUM_LIT> <EOL> rsol = <NUM_LIT> <EOL> mjup = <NUM_LIT> <EOL> rjup = <NUM_LIT> <EOL> mearth = <NUM_LIT> <EOL> mmoon = <NUM_LIT> <EOL> au = <NUM_LIT> <EOL> pc = <NUM_LIT> <EOL> class switch ( object ) : <EOL> def __init__ ( self , value ) : <EOL> self . value = value <EOL> self . fall = False <EOL> def __iter__ ( self ) : <EOL> """<STR_LIT>""" <EOL> yield self . match <EOL> raise StopIteration <EOL> def match ( self , * args ) : <EOL> """<STR_LIT>""" <EOL> if self . fall or not args : <EOL> return True <EOL> elif self . value in args : <EOL> self . fall = True <EOL> return True <EOL> else : <EOL> return False <EOL> def Ptoa ( P , m1 , m2 ) : <EOL> """<STR_LIT>""" <EOL> c = G / ( <NUM_LIT> * np . pi * np . pi ) <EOL> mu = ( m1 * msol ) + ( m2 * mjup ) <EOL> a3 = np . power ( P * yr , <NUM_LIT> ) * ( c * mu ) <EOL> return ( np . power ( a3 , <NUM_LIT:1.> / <NUM_LIT> ) / au ) <EOL> def vcirc ( m1 , m2 , a ) : <EOL> """<STR_LIT>""" <EOL> mu = G * ( ( m1 * msol ) + ( m2 * mjup ) ) <EOL> vcirc = np . power ( ( mu / ( a * au ) ) , <NUM_LIT:0.5> ) <EOL> return ( vcirc ) <EOL> def ringfunc ( taun , * args ) : <EOL> '<STR_LIT>' <EOL> ( t , f , f_err , rad_ri , re , k , dst ) = args <EOL> strip , dummy , g = exorings . ellipse_strip ( rad_ri , exorings . y_to_tau ( taun ) , re [ <NUM_LIT:0> ] , re [ <NUM_LIT:1> ] , re [ <NUM_LIT:2> ] , re [ <NUM_LIT:3> ] , k , dst ) <EOL> ring_model = interp1d ( g [ <NUM_LIT:0> ] , g [ <NUM_LIT:1> ] , kind = '<STR_LIT>' ) <EOL> ring_model_phot = ring_model ( t ) <EOL> diff = f - ring_model_phot <EOL> chisq = np . sum ( np . power ( diff / f_err , <NUM_LIT:2> ) ) <EOL> red_chisq = chisq / diff . size <EOL> return red_chisq <EOL> def calc_ring_stats ( taun , t , f , f_err , rad_ri , re , k , dst , tmin , tmax ) : <EOL> '<STR_LIT>' <EOL> strip , dummy , g = exorings . ellipse_strip ( rad_ri , exorings . y_to_tau ( taun ) , re [ <NUM_LIT:0> ] , re [ <NUM_LIT:1> ] , re [ <NUM_LIT:2> ] , re [ <NUM_LIT:3> ] , k , dst ) <EOL> mask = ( t > tmin ) * ( t < tmax ) <EOL> t_sel = t [ mask ] <EOL> f_sel = f [ mask ] <EOL> f_err_sel = f_err [ mask ] <EOL> print '<STR_LIT>' % ( t_sel . size , tmin , tmax ) <EOL> ring_model = interp1d ( g [ <NUM_LIT:0> ] , g [ <NUM_LIT:1> ] , kind = '<STR_LIT>' ) <EOL> ring_model_phot = ring_model ( t_sel ) <EOL> diff = f_sel - ring_model_phot <EOL> chisq = np . sum ( np . power ( diff / f_err_sel , <NUM_LIT:2> ) ) <EOL> dof = diff . size - taun . size - <NUM_LIT:1> <EOL> red_chisq = chisq / dof <EOL> print '<STR_LIT>' % diff . size <EOL> print '<STR_LIT>' % taun . size <EOL> print '<STR_LIT>' % dof <EOL> print '<STR_LIT>' % chisq <EOL> print '<STR_LIT>' % red_chisq <EOL> BIC = chisq + ( taun . size ) * np . log ( diff . size ) <EOL> print '<STR_LIT>' % BIC <EOL> return red_chisq <EOL> nn = <NUM_LIT:1> <EOL> def costfunc ( x , * args ) : <EOL> ( y , dt , i_deg , phi_deg ) = x <EOL> ( grad_t , grad_mag_n , t0 ) = args <EOL> global nn <EOL> ( tmp , grad_disk_fit ) = exorings . ring_grad_line ( grad_t , y , dt , i_deg , phi_deg ) <EOL> rmintime = np . arange ( np . min ( grad_t ) , np . max ( grad_t ) , <NUM_LIT> ) <EOL> ( tmp , rminline ) = exorings . ring_grad_line ( rmintime , y , dt , i_deg , phi_deg ) <EOL> rmint = rmintime [ np . argmin ( rminline ) ] <EOL> delta = grad_disk_fit - grad_mag_n <EOL> delta [ np . where ( delta < <NUM_LIT:0> ) ] = - delta [ np . where ( delta < <NUM_LIT:0> ) ] * <NUM_LIT> <EOL> dean = np . abs ( rmint - t0 ) <EOL> cost = np . sum ( delta ) + ( dean * <NUM_LIT:20> ) <EOL> nn += <NUM_LIT:1> <EOL> return cost <EOL> def ind_ring ( ring , r ) : <EOL> '<STR_LIT>' <EOL> rdiff = ring - r <EOL> return np . argmin ( np . abs ( rdiff ) ) <EOL> def ind_ring_big ( ring , r ) : <EOL> '<STR_LIT>' <EOL> rdiff = ring - r <EOL> rdiff [ ( rdiff < <NUM_LIT:0> ) ] += <NUM_LIT> <EOL> return np . argmin ( rdiff ) <EOL> ( time , flux , flux_err ) = j1407 . j1407_photom_binned ( '<STR_LIT>' , <NUM_LIT> , <NUM_LIT> ) <EOL> rings_tmin = <NUM_LIT> - <NUM_LIT> <EOL> rings_tmax = <NUM_LIT> + <NUM_LIT> <EOL> print '<STR_LIT>' % ( rings_tmin , rings_tmax ) <EOL> goodp_rings = ( time > rings_tmin ) * ( time < rings_tmax ) <EOL> good_rings_npoints = goodp_rings . size <EOL> print '<STR_LIT>' % ( good_rings_npoints ) <EOL> ( grad_time , grad_mag , grad_mag_norm ) = j1407 . j1407_gradients ( '<STR_LIT>' ) <EOL> fitsin = '<STR_LIT>' <EOL> fitsout = '<STR_LIT>' <EOL> read_in_ring_parameters = False <EOL> read_in_disk_parameters = False <EOL> tx = <NUM_LIT> <EOL> vstar = - <NUM_LIT:1> <EOL> def print_help ( ) : <EOL> print '<STR_LIT>' <EOL> try : <EOL> opts , args = getopt . getopt ( sys . argv [ <NUM_LIT:1> : ] , "<STR_LIT>" , [ "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" ] ) <EOL> except getopt . GetoptError : <EOL> print_help ( ) <EOL> sys . exit ( <NUM_LIT:2> ) <EOL> for opt , arg in opts : <EOL> if opt == '<STR_LIT>' : <EOL> print_help ( ) <EOL> sys . exit ( ) <EOL> elif opt in ( "<STR_LIT>" , "<STR_LIT>" ) : <EOL> fitsin_disk = arg <EOL> read_in_disk_parameters = True <EOL> elif opt in ( "<STR_LIT>" , "<STR_LIT>" ) : <EOL> fitsin_ring = arg <EOL> read_in_ring_parameters = True <EOL> elif opt in ( "<STR_LIT>" , "<STR_LIT>" ) : <EOL> fitsout = arg <EOL> elif opt in ( "<STR_LIT>" , "<STR_LIT>" ) : <EOL> tx = np . array ( float ( arg ) ) <EOL> print '<STR_LIT>' , tx <EOL> elif opt in ( "<STR_LIT>" , "<STR_LIT>" ) : <EOL> vstar = np . array ( float ( arg ) ) <EOL> print '<STR_LIT>' , fitsout <EOL> phi_deg = <NUM_LIT> <EOL> i_deg = <NUM_LIT> <EOL> dt = <NUM_LIT> <EOL> y = <NUM_LIT> <EOL> Rstar = <NUM_LIT> <EOL> Mstar = <NUM_LIT> <EOL> Rstar = <NUM_LIT> <EOL> Rstar = <NUM_LIT> <EOL> Mb = <NUM_LIT> <EOL> Pb = <NUM_LIT> <EOL> t_ecl = <NUM_LIT> <EOL> a = Ptoa ( Pb , Mstar , Mb ) <EOL> print '<STR_LIT>' % Mstar <EOL> print '<STR_LIT>' % Rstar <EOL> print '<STR_LIT>' % Mb <EOL> print '<STR_LIT>' % a <EOL> v = vcirc ( Mstar , Mb , a ) <EOL> if vstar > <NUM_LIT:0> : <EOL> v = vstar <EOL> print '<STR_LIT>' % v <EOL> print '<STR_LIT>' % ( v / <NUM_LIT> ) <EOL> dstar = ( Rstar * rsol * <NUM_LIT:2> / v ) / <NUM_LIT> <EOL> print '<STR_LIT>' % dstar <EOL> if read_in_ring_parameters : <EOL> print '<STR_LIT>' % fitsin_ring <EOL> ( resxx , taun_rings , rad_rings , xxxdstar ) = exorings . read_ring_fits ( fitsin_ring ) <EOL> else : <EOL> print "<STR_LIT>" <EOL> rad_rings = np . array ( [ <NUM_LIT> ] ) <EOL> taun_rings = np . array ( [ <NUM_LIT:0.0> ] ) <EOL> rad_rings = np . append ( rad_rings , ( <NUM_LIT> ) ) <EOL> taun_rings = np . append ( taun_rings , ( <NUM_LIT> ) ) <EOL> exorings . print_ring_tau ( rad_rings , exorings . y_to_tau ( taun_rings ) ) <EOL> if read_in_disk_parameters : <EOL> print '<STR_LIT>' % fitsin_disk <EOL> ( res , taun_ringsxx , rad_ringsxx , dstarxx ) = exorings . read_ring_fits ( fitsin_disk ) <EOL> else : <EOL> print '<STR_LIT>' <EOL> res = fmin ( costfunc , np . array ( [ y , dt , i_deg , phi_deg ] ) , maxiter = <NUM_LIT> , args = ( grad_time , grad_mag_norm , tx ) ) <EOL> kern = exorings . make_star_limbd ( <NUM_LIT> , <NUM_LIT> ) <EOL> ( ring_disk_fit , grad_disk_fit ) = exorings . ring_grad_line ( grad_time , res [ <NUM_LIT:0> ] , res [ <NUM_LIT:1> ] , res [ <NUM_LIT:2> ] , res [ <NUM_LIT:3> ] ) <EOL> samp_t = np . arange ( - <NUM_LIT:100> , <NUM_LIT:100> , <NUM_LIT> ) + <NUM_LIT> <EOL> ( samp_r , samp_g ) = exorings . ring_grad_line ( samp_t , res [ <NUM_LIT:0> ] , res [ <NUM_LIT:1> ] , res [ <NUM_LIT:2> ] , res [ <NUM_LIT:3> ] ) <EOL> hjd_minr = samp_t [ np . argmin ( samp_g ) ] <EOL> hjd_to_ring = interp1d ( samp_t , samp_r , kind = '<STR_LIT>' ) <EOL> ( rstart , rend ) = exorings . ring_mask_no_photometry ( kern , dstar , time , res [ <NUM_LIT:0> ] , res [ <NUM_LIT:1> ] , res [ <NUM_LIT:2> ] , res [ <NUM_LIT:3> ] ) <EOL> print '<STR_LIT>' <EOL> print '<STR_LIT>' <EOL> print '<STR_LIT>' <EOL> print '<STR_LIT>' <EOL> print '<STR_LIT>' % res [ <NUM_LIT:0> ] <EOL> print '<STR_LIT>' % res [ <NUM_LIT:1> ] <EOL> print '<STR_LIT>' % res [ <NUM_LIT:2> ] <EOL> print '<STR_LIT>' % res [ <NUM_LIT:3> ] <EOL> print '<STR_LIT>' % hjd_minr <EOL> print '<STR_LIT>' % np . min ( samp_r ) <EOL> time0 = np . abs ( time - hjd_minr ) <EOL> time0_grad = np . abs ( grad_time - hjd_minr ) <EOL> flux_color = np . chararray ( ( time . shape ) ) <EOL> flux_color [ : ] = '<STR_LIT:b>' <EOL> flux_color [ ( time > hjd_minr ) ] = '<STR_LIT:r>' <EOL> flux_col = '<STR_LIT>' <EOL> for b in flux_color . tolist ( ) : <EOL> flux_col = str . join ( '<STR_LIT>' , ( flux_col , b ) ) <EOL> def plot_folded_phot ( f ) : <EOL> '<STR_LIT>' <EOL> h1 . scatter ( time0 , flux , c = flux_col , s = <NUM_LIT:20> , edgecolors = '<STR_LIT:none>' , zorder = - <NUM_LIT:20> ) <EOL> h1 . errorbar ( time0 , flux , flux_err , zorder = - <NUM_LIT:30> , ls = '<STR_LIT:none>' ) <EOL> fig_fold = plt . figure ( figsize = ( <NUM_LIT:16> , <NUM_LIT:6> ) ) <EOL> h1 = fig_fold . add_subplot ( <NUM_LIT> ) <EOL> plot_folded_phot ( fig_fold ) <EOL> strip , dummy , g = exorings . ellipse_strip ( rad_rings , exorings . y_to_tau ( taun_rings ) , res [ <NUM_LIT:0> ] , res [ <NUM_LIT:1> ] , res [ <NUM_LIT:2> ] , res [ <NUM_LIT:3> ] , kern , dstar ) <EOL> gt_abs = np . abs ( g [ <NUM_LIT:0> ] - hjd_minr ) <EOL> g1 = g [ <NUM_LIT:1> ] <EOL> gt_ingr = gt_abs [ ( g [ <NUM_LIT:0> ] <= hjd_minr ) ] <EOL> gt_egr = gt_abs [ ( g [ <NUM_LIT:0> ] > hjd_minr ) ] <EOL> g1_ingr = g1 [ ( g [ <NUM_LIT:0> ] <= hjd_minr ) ] <EOL> g1_egr = g1 [ ( g [ <NUM_LIT:0> ] > hjd_minr ) ] <EOL> h1 . plot ( np . abs ( g [ <NUM_LIT:0> ] - hjd_minr ) , g [ <NUM_LIT:1> ] ) <EOL> h1 . plot ( gt_ingr , g1_ingr , color = '<STR_LIT>' ) <EOL> h1 . plot ( gt_egr , g1_egr , color = '<STR_LIT>' ) <EOL> h1 . plot ( np . abs ( g [ <NUM_LIT:0> ] - hjd_minr ) , g [ <NUM_LIT:2> ] , color = '<STR_LIT>' ) <EOL> h1 . set_xlabel ( '<STR_LIT>' ) <EOL> h1 . set_ylabel ( '<STR_LIT>' ) <EOL> print "<STR_LIT>" <EOL> print "<STR_LIT>" <EOL> print "<STR_LIT>" <EOL> print "<STR_LIT>" <EOL> print "<STR_LIT>" <EOL> print "<STR_LIT>" <EOL> print "<STR_LIT>" <EOL> badring = <NUM_LIT:1> <EOL> def onclick ( event ) : <EOL> global rad_rings <EOL> global taun_rings <EOL> global badring <EOL> for case in switch ( event . key ) : <EOL> newt = event . xdata <EOL> newtau = event . ydata <EOL> print '<STR_LIT>' % newtau <EOL> if newtau > <NUM_LIT:1.0> : <EOL> newtau = <NUM_LIT:1.0> <EOL> if newtau < <NUM_LIT:0.0> : <EOL> newtau = <NUM_LIT> <EOL> newr = hjd_to_ring ( newt + hjd_minr ) <EOL> if case ( '<STR_LIT:r>' ) : <EOL> exorings . plot_tilt_faceon ( rad_rings , taun_rings , res , hjd_minr , dstar ) <EOL> break <EOL> if case ( '<STR_LIT:v>' ) : <EOL> exorings . plot_tilt_faceon_vect ( rad_rings , taun_rings , res , hjd_minr , rstart , rend , dstar ) <EOL> break <EOL> if case ( '<STR_LIT:a>' ) : <EOL> print '<STR_LIT>' <EOL> print '<STR_LIT>' % ( newt , newr ) <EOL> bigr = np . append ( rad_rings , newr ) <EOL> bigtau = np . append ( taun_rings , exorings . tau_to_y ( newtau ) ) <EOL> sortedr = np . argsort ( bigr ) <EOL> rad_rings = bigr [ sortedr ] <EOL> taun_rings = bigtau [ sortedr ] <EOL> break <EOL> if case ( '<STR_LIT:d>' ) : <EOL> if rad_rings . size > <NUM_LIT:1> : <EOL> rsel_idx = ind_ring ( rad_rings , newr ) <EOL> rad_rings = np . delete ( rad_rings , rsel_idx ) <EOL> taun_rings = np . delete ( taun_rings , rsel_idx ) <EOL> break <EOL> if case ( '<STR_LIT:m>' ) : <EOL> rsel_idx = ind_ring ( rad_rings , newr ) <EOL> rad_rings [ rsel_idx ] = newr <EOL> taun_rings [ rsel_idx ] = exorings . tau_to_y ( newtau ) <EOL> break <EOL> if case ( '<STR_LIT:b>' ) : <EOL> badring *= - <NUM_LIT:1> <EOL> print badring <EOL> break <EOL> if case ( '<STR_LIT:o>' ) : <EOL> taun_rings = fmin ( ringfunc , taun_rings , maxiter = <NUM_LIT:1000> , args = ( time , flux , flux_err , rad_rings , res , kern , dstar ) ) <EOL> break <EOL> if case ( '<STR_LIT:q>' ) : <EOL> print '<STR_LIT>' <EOL> break <EOL> if case ( ) : <EOL> print '<STR_LIT>' <EOL> calc_ring_stats ( taun_rings , time , flux , flux_err , rad_rings , res , kern , dstar , rings_tmin , rings_tmax ) <EOL> strip , dummy , g = exorings . ellipse_strip ( rad_rings , exorings . y_to_tau ( taun_rings ) , res [ <NUM_LIT:0> ] , res [ <NUM_LIT:1> ] , res [ <NUM_LIT:2> ] , res [ <NUM_LIT:3> ] , kern , dstar ) <EOL> gt_abs = np . abs ( g [ <NUM_LIT:0> ] - hjd_minr ) <EOL> g1 = g [ <NUM_LIT:1> ] <EOL> gt_ingr = gt_abs [ ( g [ <NUM_LIT:0> ] <= hjd_minr ) ] <EOL> gt_egr = gt_abs [ ( g [ <NUM_LIT:0> ] > hjd_minr ) ] <EOL> g1_ingr = g1 [ ( g [ <NUM_LIT:0> ] <= hjd_minr ) ] <EOL> g1_egr = g1 [ ( g [ <NUM_LIT:0> ] > hjd_minr ) ] <EOL> x1 , x2 = h1 . get_xlim ( ) <EOL> y1 , y2 = h1 . get_ylim ( ) <EOL> h1 . cla ( ) <EOL> plot_folded_phot ( fig_fold ) <EOL> h1 . plot ( event . xdata , event . ydata , color = '<STR_LIT>' ) <EOL> h1 . plot ( np . abs ( g [ <NUM_LIT:0> ] - hjd_minr ) , g [ <NUM_LIT:1> ] ) <EOL> h1 . plot ( gt_ingr , g1_ingr , color = '<STR_LIT>' ) <EOL> h1 . plot ( gt_egr , g1_egr , color = '<STR_LIT>' ) <EOL> h1 . plot ( np . abs ( g [ <NUM_LIT:0> ] - hjd_minr ) , g [ <NUM_LIT:2> ] , color = '<STR_LIT>' ) <EOL> h1 . set_xlabel ( '<STR_LIT>' ) <EOL> h1 . set_ylabel ( '<STR_LIT>' ) <EOL> h1 . set_xlim ( [ x1 , x2 ] ) <EOL> h1 . set_ylim ( [ y1 , y2 ] ) <EOL> plt . draw ( ) <EOL> if os . path . isfile ( '<STR_LIT>' ) : <EOL> os . remove ( '<STR_LIT>' ) <EOL> exorings . write_ring_fits ( '<STR_LIT>' , res , taun_rings , rad_rings , dstar ) <EOL> cid = fig_fold . canvas . mpl_connect ( '<STR_LIT>' , onclick ) <EOL> plt . show ( ) <EOL> raw_input ( '<STR_LIT>' ) <EOL> exorings . write_ring_fits ( fitsout , res , taun_rings , rad_rings , dstar ) </s>
