From a1dd588679d5b058ed9ac2cf303af8e453642d07 Mon Sep 17 00:00:00 2001
From: Irina Ene <irina.ene@berkeley.edu>
Date: Tue, 6 Dec 2016 20:09:30 -0800
Subject: [PATCH 1/4] add benchmark file for passive models

---
 fsps/benchmark_iso_color.hdf5 | Bin 0 -> 44336 bytes
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 fsps/benchmark_iso_color.hdf5

diff --git a/fsps/benchmark_iso_color.hdf5 b/fsps/benchmark_iso_color.hdf5
new file mode 100644
index 0000000000000000000000000000000000000000..192b7d06f44601f65b329b4b049603c0c707a17b
GIT binary patch
literal 44336
zcmeHw30zFw|Nm{D_L<gck*Soe7)eDQ_hW4l4Jl=*$x{lI6bWH$iN_KlB1TBcl0Ey_
zvzAaHj4jy{MWX(nY3}t?|KG#&d>{3Dem`?xujbry&*z@cx$m=d&$*xX`7Chl(@#^)
zSdEYz3<jY{FzZk0pE~N`?vP%g<8(ToI+3Ho1S-ssm!8NG3__B|rmi!o{5cBJ3pvz{
zLk9QlOQ0m3KOIuH5yc$oS;NxNsj~3@Is&eJd-5dD=u15)hYEA4c1A^sUV&ai+%Uom
z3>qI3NDu?uyoN~9==SME3;xBBMYY$)?OKBX>?@Ld8H67G1qF`|oD?>ZFhWI=HY)j(
z>SWYwqV+3bEOL7NrTTOY!j;NLtv;X8W5R+1A_BvMr6NDfPwSifBv3TYZ{nEn;OL2g
z;gNqaKjUxmlXTJl%1^Vtd^lWm3GfMyiXJmQLXzc=o>%#s{3Lia&M&k}<dnaY-?#Wl
z&ZS7#pEUT9Bh;k$Q^3nnBaob&^g<U(nhq5(U6R7$kXW0N0|<L{!oH3gr<AC1OkQ5T
zJ_80NuOP>#hPDh!Q~k|^4wZ}KB1hZB_rKbUay^cdBF0RN#<U;dC&7#_Ef688YkK@#
z?SK(HLZt!z(j)BG3k^=>@aZIFp{3|{&6HFS9(+!HK>TPN&S@o;leBNy@#`TVM<_LD
zM_=2cy%dYwMv%Dl_VlRl7*vC)KSvp9-zA9Ph>1~CFisESk~S_mMg>Rw9vU7PHLgQw
z&}i(WR5v6f<CEmD<92ee?PO=mbr{mAQ@2h|-MTn*uy?TU>_ByBNxzXCRP!b}b#Us?
z$??DR()#T>ICS8CBS+sC9vC9&o`j0z2M=9oK=d0EHW5cZj=PPkO@A9V8y}lq^(j;r
ztsl!J7wG#{sr%^u(#v*KU_=lHlW_o;6cj9d@gn}e`_voV=i2B#4Jy5Hc{Ll|=ho;x
z8sEnEX*ashz0rMi`5WKIYII-!M)&FbNjYge>3n?}-KR(0M?a6p_z+IdkYG-mDN)gF
zF-@OJH>E;4ef($+Z+zsWu)vA<4~jD;f`jsuNPQ4XXZY<!*|`h>83Hl{WC+L*kRc#L
zK!$(}0T}`^1Y`*OV+2GHwlU)PTT-woK5j3wf)sylHE>=~1u486v|h*Y9VuGs*UwPz
zJt=k^E1DZ!NfP-%lkEzsNbv+Ghho!@BwS9j7N=E{{BGxur6hhL#YIOx&YV|6iqhYH
z8F1?}$)5_v-Q8+QQQT>Zz}t0XWK~-9nOy|%N>@*LK1>ccRtYw~dh#H;zNYm(A9)aV
z&E`*7D-VgksTKA%QUK0TUEd|C3cx?yMvco+1Y&Yg%=`I@Al@`*`n+~Z5XpEJ(_5$n
zyjLz$^<OE0*ved8cvcA#CyMv>{G<q?=YqRg@fc2>8NJ_K8M3<%I(IT$1&9%H1I<YW
zWM>&Gh!-&+f2>hne`i$?d?|L1{-6p_y_dHoTMY!8ta_XppbkXrj0uO{ssqpIy<kMN
z24uJ3nH1P)0z|ZYtt{3A{*2_4_6xK?pf}~rxYpXhTi|o-{0VK~FC1Q&+>QyMI(~+$
z=Q4pe<#fo5JSLQ{9AMf;g9XIZaffn8vVcErn2F;~7UXp6lccJp1454nZ?{D1fXI4%
zi|faAa50;^t7iNoQWVZwuO9Y@BqmzdnS?$f1#!F~9lJj!3AckKk@`<aFqF40y7!cX
zq3z@Jwmc^Z_dCzz$ycOs$c}_6n{rZc%ygEn##@qeHnZGUy@C`aPCR*E@Q&mM8t&Q7
zt|SH3rktI(sz|{hzYQDaeIz+m6&F08WBrD{^luqaL-M01xfh<SA;ktQMzt~gOmgnm
zY#Xlpg(Qyo%`>b0Lc-We*K03pNrI!1+ka~vDHi&X+Hts8tH?gm{v83_fQe%xddWfL
z)LQvBOgZ4W=63P3AwZZ}EsWbyOA4pPWnS%!?Sh$fv2(CIMAlkWXgDhX*Jb$F-b)mK
zdG=Vp_dOJWvn9vD`ZLz6XqL;>eM-0(Ty#{2rwn|ZK4Cepl|h_to;M^?1^DL@-1fC*
zfY?8~%hDnShz47_3FfIn<oX9&i`uCHJbknE=|we2<S3@U4pj&K+&$HWtuz3v>pWdV
z8o*^ORq^uC1hGS`s+XP?L>~Y2)G$vAM7*Sg&pWjN7TNogy_w)NBB|unN+t+~*Lm}f
z-6sVzc1@a^exD>F7<${6Vz^j%w`a^fl9*k$p=;n>62|qMJCyf;BoxDMq`Y`ULPUFE
z-;t%H=+!Rfi-FHbLcVI;_I=Muf!|9x_Pm#*U=1hlP}Cce$a9{!cko+MI5&6Ta{YIt
zpv0%eCDD6Qka;ltd`qmiWgheItouNU(iBs!Uim-@)~8txsQW+?QHxA^JorfR2WG9f
zbfKCAtrvZlWqcw9Ck$^~^{pWVx3?|QJocF+Mn+w6+W&<VOjsIP5Qg>G=b?EEJM6y(
zYdhtpVt<aB_~Gytly4WKxKLP2^81A-d#7W6-u3xr?Lw^2XNJS$%H=><;uNT{LmqfJ
zE{W{P3Lx11M^E=|ijX+!yv63bioo4jd);!55^!e}40tE649rtA_ntbc3<BkyV}w|5
zBEC+=P^>rRf*s2=ZZm*CTmDJGEL9M>O)7rZRt=bninC42)qvO5^MUCqbr3z#3cu4s
z0|Zx<bn>+{LFB01a$2z_Brdak{p6$;aK%ejFG$k{k)rk55pA)a;^k`tVyj7k-Go<f
zmVY9llhuYDB782_-2v+Q`21m+=hEiYk`QNJs{Epk6fGYy>Rl!Q{7i-HCBw1(PIn~y
zta02*-l=bUN*=g9T|^Q6v7NqI_h7DA0m|2Bs(T3(0bUGSGusBoqtjRK>Xaw}v%oNU
zZ-O%9hdOn**INa+fld!MomT<D3;#g}z3@4QmYAKbVnDfxvkP}CjzhO%Q}>NjgZy@Z
zuisRwf$-p*0)utxkgxn|))kfpa0ZSV%^QmCb##^Wo30v=IA@u<%~EyX?oOEFeO?Vj
zn@5`!I%`6H<(`lsT3V2pl5_ChX~e_oK($r57Kl|{uRl%FhK%O{rPUTVKIJQ(E{(%b
zS!Jv6BH|>@T(Y4B3mBKr7>57Og6x%Vb!<<ufH}E5Z<;BNPwk>#?g-a`?9W5O(hp+X
zZCO@NIbC3@R+dG%=t90?)ZoHYU0_};;`e^43uQ*OcV8Rmfp2Axu-cJ&P<)f6p}Rp3
zxWzhgryuG8-(&pqqn-4D<7}xlVwyhu<9~3n48Ds%<L|cg_a_0T;cY=&X?!9T(%%W7
zq1a!Jc@6aJBZprN=ymxrYQdk*m)YpM*Vl5zNLA2s8B{zsR(hRYuNL)aC|5?s89o1U
z5&u^R{^)xkVJ&@(HvcQ0Mpnq5jX>l6?L;Zyw`!=#or=0rVSg%=*lI`)A89~HzoyZ)
z69mDbt~U;({agNP|K?I1o$j|B3+Wv|m4oiLVyb|&9HXmLu8nhJ4ZlWQilZ_?4ZB1_
z=kGr^jiYPp@5^zh<WR}Sy6+tH_*_hXSo{-NdVKc%!9YVN$F}>EjGE|j^Vyp)s88Tc
zZ*L;`Hcfxu&!xVrhknJ6u92@78k{sA9@!>R#Z8Zgr0vsq5NJyvIsSck__mOitm*NP
z*dWk&a2=#t|9yBQv~MmvB(@GT9-)7;eZXcX&5M@_D4!TtZBU2vTcX3-@KHFgeG%j`
zo`LHSNAE2xI)dS<SySJ7DS<#~@UxY#ab2Ql)3LI6W#HAgckSq?0>aWd^Nu**6BqAH
z>e3z8oi?=C`eh%6>5&cx>~US<S=&;_qpFa<zk|VBo*HoHn{SBypa%T#b*|I1)PXbp
zQ%PMv4G_Ct3-9t;1Ip))Kk1UH3FRACOxEznbqTlfnQPl?L%xqmd9R1sAT-GHcK2X{
z>)Yv?vyz#>4zb(BDq`aGEN!J4CM2p`gloTLf-fvd4!FmJ+*x;D9BQotJe9hHZbBW%
zjS3yCr-AdfG~0m%4mgkN?5VwTCB}F4=Y06A3!<6Bmrwg$GM{Td{Pkfy;JNErw_xi-
z>~q7>a$_*Oudkz(r;qdW^q?9$1K>n`nWHer0E)L1Zd`TV0LssFE?vYl1kwC?h7Mka
zkXWnVbxpb<@WakWrM)qPj3u+RXREV;@3tksLvJ>Sht7z7Jeds%*JiGXPQ~kQ#}2x3
z9j~V?GJ0!e1e{s@cNkuh1LBXl=iW?Hfb37#Mm}Gu1fmI>gIw%!-aS?)>*76JXE>6w
z&~cm!aE$9>vZ_^p=fl`sJ&ysL$cief_Nu@);07GqkMr-zjOF)QsX=z!ABlZdVm){-
zzmcwv>wz_2a(5@FgD_yRMjy5Yh~(MIS*LM6{;pL`&~!}@$}Mlvqdm^IBTm>CKh*-}
zgwLL1GI5^Wy=>Xro=nK!(?TyU2kT?(rnvZ4xc*mAux$Hs7L<7o>giX>f=JVJd+p&m
z;49weSbPBI*HPy07Jt$KvFW(wHl!{HcfFjF8?6Vfxji{Ciuw@gIYlgr*9T#~Tjy<8
z5ck<9t(SrU_-3sh`s|ni@Qis8?VSuE_O{NM)9Vd^oABz%K~*+5U1`P04Q4}m@b4$b
z-N5+p+TX_XGJ^boRy)@vA)a9-XVY#OLHRk2S>;y7z=${4HDRnVBviUzy&^J({1=1z
z9?~|!=axSo-ogYJ*CVwZ`kP=sxs!4<+5|%F4tri6Zvwagka}%`3FHrX=-2tE2}t|A
z8Fa2c%wTl=VFoE{B!?MXxhy%%D7luDTuTV4mupGMwWQ=)QgR&Wbt4WX&yl8kaH#t@
zQaJ%7w}z72OUYfP<Vq<y!ko&_9PI1MX-?&5PUUA#<wwb_q2%^ba+fK&Qc8}npz^bz
z%4b34XF)x`1(lx#)WfiblG{tkU8dwpDLKMYDraO#<!4FdXG!H}N##e$t)b-hQgW9m
zxl&4wXhr4MipsB*^nOBi$PoCo5oo^mWv;Ch>85`_WS}DYeTe@z=MS`fdgFQg&#<GH
zEdBi#4?)Vi2pb;}@^^R_68lx!E?=b{P<m!uANl(8V?A<+kpw4t%|n73eJ;ID>rt&w
zX?i_U+vBF!BQ?csdOcFx=BC#pwS{haJ(9L3Z5j0^4SwVZ<$C*k>5rrxN$2T5pZ8MS
zL8|!s;|M{-HYi_nNB)ejP}W^ql3c1lKU(J&QtN}#-r3wwA=CFijX>jdZn~GdHd;4E
zjzYAgtn~V_>pv{GbbVIBGL#<q)Vjhy%=VYNoqnHayx#wxzE4R0chUcsC0vys|DTtR
z*phZ5Z5z@WrQ3;w$3o{P<+ji*3=_UNn(sX^+gqyk&&T)lf5>ah?km;sEw51!A-xv*
z^_hRP&ivo=8XJF){n7XieWVKid(^)cMt1Ul6oJO$JKgHq^lbRRKif{>&)!d5gQR+M
zDfK^Uzu6<Df^8ZWR~Gw45NO<fov7z7>)jB5rv45up!y5_e=UQG=ThaM$JYdE?5E{&
zBc<gmi)vWjUnTZu<3PqlDFXC3&^&x<JET-|JJ{QCojP-+`svT}l-ZPwPfhW+Wqj&i
zrH9I%Q-**H0T}`^1Y`)v5Rf4tLqLXr3;`Je&5MAjn%{Sk?kiFt)KFTutb!DDwR!p8
z5Bby2riZ0lRgthwF2}P9d6}BT?}}!hNaDTx2%WJtBqzvO?dYb@q>#PI`g(ulTZ$)q
z&hV%s1zRro_q8Md9EPktzm-5c!;{s7wj7iT)~{nFAx~C0w@VNUd6r+kh>tClhirk*
zSr0Sh$-4H<3$sxGqO<1^KPHBz>`58r3LqrU`o24_2t4nWCw7J)Z`ygl&5H`kApC9E
ziWLjcPIX|*OV2(iLuCGGEiUpfxf_f747h@Jl0N5eBpES4u+L%+*BisvVXrryLA%xa
zJ4URqRE6>djAL7T&~CEwQljgA1`zx2IjYV@JJ-E#JI!9J0`s(8(zE->r}ivdR-A(T
zMB@}ElB)r{=-48cIOJ_6yqxF$7J1coj?WmEf_%k4W=}Y1gZ7gK=Pz^-BVW<EXz1?A
z+Q2z6eccf&CW!Vw+_Yg9@~GL1Hp{(cg7{L5#X}_)eB9G_tCKF;Q)b@To#u;SaptqN
z7g)fyH*-J9&;hvTai&9<4hZsG*B;qZNQypWq-NP5fBM#veEXm>Qn0~z7jyGV5_+U<
z_MBW!3J*;lwB$<#$ua-9L=cblczNZj8x@tLh#lPD;^qgElQBKN<lRS-5d8k&OD*zA
zJ9V$MWFU`pxa!b$=`|!g&FXu0HS%W1DA;VD`-Kz|DtoSm*OHv>r8Bk#*O7eIu-l!V
zVcao1c}zL-PCXXpkR^B>4ElH!5I{_L5pfv#rW{3ug{Gx)!1u|?^!;5PKsWS+=P`K@
zD|nX_v{V35*BBQif8@!}*<Ta23hR5r`qo<xA`kVb{?=McMG%5(M-Al7^HeMfRnrwA
zzbdJAY6aRaYup^*y%*chPM!O1Wys6DX&i0%0_~OyD&jWn!S>K;*opStRUmT6oo&{!
zXuo`Gp;g!m42PVmb)3%t;Q;?dZ8%uJo~&yFZ%XWzH&*$_ATM{+dCL_O)d9R7|EBAL
zcCj%&DciWn%Z<-CC5q7mf#Kx&-s_M*AD9@Z+zspZz*wbSw%R~c78Uk6rVV135qj)F
zOd#}(4uo&LL=tVEmF;+lyy|Yb^E~_?lkhA#W>tO}35*X#9v-hqV#~C_nSGJhe&bZO
zo*4HpkdNKlfB1V63MxZ~^sgd0$8YZAMSUQ}{E{(xvp<sJUhyg87FClXg^7FD+(v%8
z<&C}fcOehhKU(WrLJdizZcSS3hy5{r*pnF@kl((%bjtd0<VUZMJ2Tc4>viF_g@Yv5
z?Jp;f5+iTDPrTuWeRX7cpL5zvej`AzcJhq-DadyZiRoBYN}yfwBKtLtXgA(r+n9`K
zIe-h7ACzoB9&yE71N)0|K=6;XFz$o<AFT9oTa5LQ?X7l1<u=9(*yHm^jNACk46>2~
z@mk)>gq~Q>`$Ai79H9vO-ZQ@3vqSslK8rMW;JybU<4@$_&$tglYCa3*Dnq&6ki+IS
zDv;gpOv$7?6(G8`3T>N+JmWnlpQl@>LSp8K-FuJ+FPxt2^M|E6l<OMnv&=PsUw2}~
z<1h@vKcsI;LVI&TzmenYu-?4RrtL?&X681xc|}(++_RNE=``}HM~z?I?H-PA+X-T#
z9rDVzjngq)i2DRQPnvJ|^fM{Q^S)i@fqd$R34_OeK|VDw+A<d-pE{iwKh#7Hm`~sP
zJT8(0(d6!h3cmQf#FSgt?_m4B^x}<ExB~Ew=g9T=gzfy${;2RJiVzvFX!@wuO28jS
z#9JR!0*+YIG&e{Y62m{m&b)+n-i00RymMB8d>4m=y_(ogKLpqL<G94VmuEC;s48Ul
z%<1ch`#kWN87=)6se!o8U*Fge?Ydi(rnlRn4v7~ZJ#X6y$EUL<)B0_}_L~vX=?Dw$
z>@6x!1h2yV0Pa|unrk9|`pUymJEL*mhfD86=BQzO;f7sD4G^bsAtf%Cv?27$(2)~H
zFd?Voe8*iUnP~Suuz#!$3t)xys9aALB<fy!`Xm+k)JbI<t3I+oyuU3wl&b?7t9C4%
zJp;p*`cscw(gCM8{hr@!qYLoS|0ZvyE@YHUyTAKB@~dyfe;QSx3*zU4%-_{vdfjqu
zk+U8!P8@e0v{4V3w=4&byN>+o1>Hw_n9KOoP5m~|n6Kyb`}Z}{-&;QVe_xJ6C9?y+
zFUO#gvw!=(9R2;2!T-J-{k^iysPD`D?9}G7qwpOnEql}3QK+FJ+TJ!hT&gwnD{yp;
ze7(@%r1|h5{G^JT9uG;|r|~HJ+jb?*wSRcVu;#+!TfQFqUyDa}aC71DEnkl}UaIxq
z*I#@?nhOuf*iQEsCRGc66%X3JuJQN&|J1%t^1Ye15Bpj6`gDFxZLd%3`JTPLv`+v2
zFK4%%x%Kf{Qm`a_(3v$jZ$CfAZ!8Pn_m9jk4o$@OcY}}H<{RREK~wE(&*S`!tD5Lt
z(@O!w9(q<``8aQXnY(uHT-@Jm%n|DXecUf-!MlDJ_9#JlR6ywTK-|BqTzBeyEu4Q1
zPPQpNjN!Q{ZSO@gAkyltEyD=;ik8ca`aZyYl3e>6CC|Y5w!3b7Q(JY2<T~bVd5!P)
z#~v3RTPfMs%EIin8@|uB(5amF1?Srf!_2qt)&l;NG+t$hHsoA?GRt3&3EA84X0D55
z0-T$8nQ;Nb%tXyJ<SPoFWX~Hkmjy&xT>v~{0jEn3ui?FrpBS}e!x5eil=pt!+r|p{
zi=A&qB_!xVLe;z6-dcL_k@tRB;cPwNl*MK_ywrm>;=IxV7k%&*U9gz73*S#u&h<HJ
zVgODv6zbacHGup@Yr+@JGXVHQczw}T0|@2ROvo}Z1YS>PCHF~&;M=9_(C$J*V6*iX
z=$075N9U^P{x1z7Cv(xjxz=nb`%o7@GM)|0gRR%M-h;f+BjLX8Z`lAL%45A(;5_@q
ztMe0XV?BwxUp?!I^X$U*(`_!x0oR<{Mr#n(&pEwrdF41iE?@rq^#ZIX#Uo);JK{X~
zWm|8u62polfBP*s|F|=0&t<G1BBCntDbBYdXFeX4S)~Gyr`%`TA)J5g6)L@n$NlA&
zuty5}BCnS}SE$ax`mh~pzq?o+n73mugsj&9W^Am^^>E~m=ENUOw9|t82S+VD6tyAI
z_H6C?BINV-DpEe+hx?5s$8FEO!bBcsyXP)0Ebvtu%_};|g4JQwTSt;QP}cSFiHE0k
zfUDL&gY?h^W~;6zh85}pk;*S`&)0+8M=LBh?9v0V<-Qxci}WD$@?p(rQv=}IhW8wN
zzyL(fWr`ZT5a(471=(i}AzMKy^~`T<$em{}MB^|Ud|xJ2ueCP<?#lg#!~Kyb+Aip~
z5o?XWNxWr8)MF#a={3dZb30>zZT43C<Bfs4IAk~boiXsEVkU>bFa|-V3quPljUiE~
z^OUDnCJ<S<W{JjR<c*30N0z3T!0M|9I={If<B3Z7u`-_M-yex(cgPTEeguBZGx#x|
zkbXaDYMvjhr>S{<w4SEs`O$isn&(IB(XDS@^a@VHw$^++KVpnj^Y`cL|I<9bAKAzA
zqoq$Fnd8WgUk`!C_VIMDVo>dZMx?3vdpXoP1FcUM{xt;X^@YYff&X-ULGnL|9%p|R
ze~-@Z-@)IDTr5TQ&&TZw=0EdGNc*vr*uc~4MvyMm`d8PD{>(37<9%0uG`<sYQoVG)
zmd0eizX}44$9MW`M_Ft`1e%M#$5<!Tn@g$xQTvTtDHUwfu(-0=FM>eh_UlAFcgD1a
z05mm!52i@<$x&)Usdxfa4tjibo!(F`AwwF^$fR!k#ezTM?}ehhisYcjfo9_GIdro}
z&!&!DIy-Tt`epn*DGoA!s4{=3(v~3ml_4NQK!$(}0T}`^1Y`)v5Rf4tLqLW=(;^^x
za#tT>UXY^wQ(ijGM!$<m!wO8>z9q#AjMX(SRG>dZ&Ti|7ccgHBR%@S|?@2z_7h2y#
z9$2@}fd_|GBkyQlne8X!6Ip0>Rb2gv6rMb%(Bau<5|ZtUP9)Zn+&y=$B@ZSbpEZzs
zvJUxO*Gu0lJ%oM}^Bp!kaz?(-=RZ=`cSio!O4h??H_-p-Vq<s|s|eZM7M)79MSqEY
zEw6@sLY_}4lM|nU_GVhsjwB;rPgoG{H6cj_5*L_e9a2U=22ZxLw(miImy2C}I}bwt
zmlfldY&(Z`XlD8oQ_;^Ub91mo$H!<dsBafgl8tt0>o0a_wG-_K{rYKCS~CDtY_sxb
ztD`^D?5pyV(9gvFx%VeO(FAVwrz_EF=wITJhuL6VCVb2uzGWHu3FO>Anw%Mdd=|rz
zW>5B_-^HhktLF?wURTHmTXVEuD|Tq>IMxUK3NC!BY;u<cvFoPoSUnZ}qUPsbZ+Dpq
zt20YaX^cjDwZWqNV*zMaDDRr>wHo~#d<-8fUxoa=2SeUgzt;m{@kP%KuXKR89}#%m
zFCYCknr%%BDj*5%CkDY|?vtV&+Hj*I`bT}{`J~Nx<n?u2Ft^R)rzCMg>EScZYm#%G
z|7^wq^n?1+Af=<%TT)c>IZgW!^2RcXJkD_6lf<XoLJL+E$qyMjI;0=^KW*3J;&}9z
zDX3|&;P)vtq;Mm5l3CnmQe-sa&Vuc=q<Gy7;*bi~XUWP6`D_fmOY+r2<REdx;E5eP
z(LZETN><h@<m)*TNo(yD0BjwvFfSsXX;{jk9JC|k$Ln0X@Dlx*wuoDD+ZO#r+Qr3b
z?m@oZ;`uru21)>$lhcdNpg)mG+6F!Sk=OOAOLA>1v=1GBJNqvB9}(@cx%e~(>p$S^
znScf8M>2ebUD{`?zjWsVTbOA7c~MbiF4`F~)9P%;FG4?&Emd1gxQO<J!+phk5w;WE
zy72qq$lF^d8a%)Pd0<{$w!J!qek}~G1}k3FLO;ECV^6n69+>AH)`1Jihw6E2zl$o`
zfd*(i-rZFTLN^93zKJ|b@!QQ^4b!zDR~&k0)HSsGT)-hxpCEsaJ#K-s!4Z<!R^#l>
zK1B+4&WQ_*DI^KYMVZW;yCk9dTicFFVp5>Eb+G6P`nOdE(?U)3TdTJ|?9Sv@r09B@
zmt4kc^gFzSU7h@f6pWlVx?MH)$7x^M)I3AKwL^7JTbUs5Y>)hEW%PS0d~<m5Vt?#^
zk36rB_d<Wf&wab6E<+xlVc791<X6JV%PgXE9Z4+fcW*{-^yl&}N8y4k+UFLw)$u%q
z_POpW15ds~er?xK_03b{fp2BDpyN$>AYK;4wAicwk<K=6BDZ6GvMaiqdmw*uW|ZF!
zA=-s{P91169P9C-`<-(Sk!L&5cA)kFl*`W^H|!1e@8D~W-kxYDIyNfVO&k54^4h=e
zwhH;N1%rH7Z$&=s^tr7S_NfAA@WRv6ZP3q4aZXTPum<`&cHQ~{d6)3M<?T0}&|h)P
z7Jk`#tlux>qjsy1x0dIXS}x&H`WZejWuZOkjgJlkMR<Q+TEu%TZD6WAz1YDU`(vW+
z{h801z^^sF-4^|k#Ii(nRp?(>{K&4Y<8t)#Ii)f_WYPzcpHSk`(+d5>zN%e2AqDx#
zC+@BF=HU3nbv&WIA9;KiCiI?C^@ZdFb`R*|kK<Rbgi&k-<S|CNxaIjHe{W{OrSWJ#
zD~=qRLX>0scYiqkCE}L-=FwBr<=C$0<~c8KgZ$E>Hy(B(4D+{6=<kYr<Gy`+4=_TW
zqSc~TUc1pBOF`sgPhV`G8&9`-&%$x4+l8f@&R{!z^n0h;DD>AG<g-LxNgc9#h1u@h
zhkRnDjsAXjw4WW2*3aoFw$sZ+Pac>d56`{N`cDE3dn6=oYKuI*^}XgV#C9#5RTDY#
zJ%*EZ2oiqN0)o*ki)o|<8KJZ8*rcO<?IiQmvYsq(nw7XQqa6L7j@D1LUxqxqyt-R2
zJD?rysgg?$VqI9BV(Z|NDDit!WLW$V?V^>BPoHLk{&VhKh`2NrdB*nVimrR3{jF`{
z#^YkNcea=B|2|w7xbnmEHXK8~UO@PYRo?pGG`p`$tfL+zW?hwg5s7xsMrpw|S?CWn
zYC}QZc0&NIF4p^I%lLau``@1+8uRx)rhH!`{k<h;)%WE%RC4j^@5?c$<j6_im!rR*
z<|cn%t__u3w&weCKRdPg_-gE>(y}+b9mTi&y~r(6t)X9mqif{rg$5_hhleOps<`R#
zkhFc;t|a4c+m$pIe=l@-bK&tVf3Nso+dd*wn+uO``FqT*Qmy~K{=(kSTzE*vcDldh
zQ}!o+6%X3JuJQN&|J1%t^1Ye15BpjCJvzUp=I_yZnu5Pa*v}}f&LV(0t(}GcZ^#Gp
zh)?x(!}oXxUSPE;&fiXI4C^`{=jVUK4mf9n^DyW0F1h`2K0aX@m&ixH-s?*N^Bk3d
zQ&D4h@|ZG+w%9UOr{O*n{qDBcv}XWQ?zoG~O9l|bxf>it<9<0?c>6cuJX~<#;>>N4
z=wEPpSht84_#Uqxx+!KazTc*9(X$ze`{E4!U{!$oPY8N-%B;!J1ir}?-(%sp|Auac
ze2FRUm*Q9Z`Xuthh$rk_m2OOk^eQPx$Uy&&A6CAPML&-)!ZYW!8w+w%vdeew$N6{W
z#xr25gY!?P&zr|!*vI|K=1?8<Q<&lsyAR_l-Z{+gr3YoR*+;T#kuP?R6MODBzJES_
zK6CqI17Hqfg_^WO-rdfS#eJTmpGIH)+{$z|d_2Bq#8|OD<lh-}rj-@?cQgo7t%x>&
zSUZJ?DZ>nb$titvqa5XezwB-mYzV>!7dj5zWeCNqBF!A#jKEhrXNtWN^7hQ0oS)Fs
z7(zw2LhSi$;Arid=P<?u5^CmqJTWx_rbeF>Vh_$YZk;evlg!hyD!&vi#{F(C-w)|L
z3HQ0_+fVZ~3-<@{a8V2>P=Z9Sg}QR5abJ)Yr{9FR;5<3`UVLV*3NZb1_F9a^b&p+c
zd2LN`e$hTT>B0fzPky}Pe##r?!$)nL*5%`TW9tE%*F)7Izs7dR<ag@8i!(`^l!Eh&
zez)d!22IFb(|-FOA2fm37V~yfwic9c8?`YxR2#D85`w%8aUR^AxnRt6+_%U-#b`k(
z?pJif=XTLhoHv)~1y9ATT3~C%lf%7nKK88R-ib$aAawq|g4$lX&?akT5%a1pgnqL8
zvY4j_`FRzVThnpAvD_>-*Ige36Y~0LA4Oi|(ivS3tVF+=bIM!ok6?qqHYxQs7y;Yv
z`Q47{$hXt$`sC(aV_<KL`0SWz0>w+T!}HPKBRhyFIAMf5zHLX-21g)|uW+o7!zg25
zw&bbzQ!)kKypP_Bk;brE`S7i})9CN<sBq+pKvT$9-*5R`$prY(=c5ALTR=|CN>}fJ
zCh&2C!eFnFW{@53_r&bE2}s9QnV-jcf0iCHKabM(Ap4ae@b5(6XYlvv_migP@6mcR
z>!GAqA{yfQGoMGhne;R@e~-rP-@)G#38gjm{rUR8oxjJ)l0JoGjw3sMJp>xt$J3)x
uZX<j1rsnU(QtJ%AUJd-SC#KgI8uJAH)Aa?(|0H^x{aO4yI=}Dn_x>MF!b*?;

literal 0
HcmV?d00001


From 1d52f7b4b5d67cfe8f1ee044aac0f20cac6dc582 Mon Sep 17 00:00:00 2001
From: Irina Ene <irina.ene@berkeley.edu>
Date: Tue, 6 Dec 2016 20:10:19 -0800
Subject: [PATCH 2/4] add unit test of color vs. time for different isochrone
 models

---
 fsps/test_sps_iso.py | 267 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 267 insertions(+)
 create mode 100644 fsps/test_sps_iso.py

diff --git a/fsps/test_sps_iso.py b/fsps/test_sps_iso.py
new file mode 100644
index 00000000..fc36dbde
--- /dev/null
+++ b/fsps/test_sps_iso.py
@@ -0,0 +1,267 @@
+# -*- coding: utf-8 -*-
+
+from __future__ import division, print_function
+
+import json
+import os
+import numpy as np
+from numpy.testing import assert_allclose
+from .fsps import StellarPopulation
+from .__init__ import githashes
+try:
+    import h5py
+except(ImportError):
+    pass
+try:
+    import matplotlib.pyplot as pl
+except(ImportError):
+    pass
+
+# As far as we can tell, Conroy & Gunn (2010) Fig.2 uses Vega magnitudes
+pop = StellarPopulation(compute_vega_mags=True, zcontinuous=1)
+default_params = dict([(k, pop.params[k]) for k in pop.params.all_params])
+libs = pop.libraries
+
+def _reset_default_params():
+    for k in pop.params.all_params:
+        pop.params[k] = default_params[k]
+
+
+def compare_reference_iso(ref_in, current):
+    """This method compares the values in a previously generated HDF5 file to
+    the current values, and reports the maximum change in color between the
+    input and the current values.
+    We assume that the reference file contains J,H,K,I,V,B magnitudes computed
+    at the same log(ages) as the ones used by the current model.
+    :param ref_in:
+        h5py.File object with keys matching those in the given dictionary of values.
+    :param current:
+        Dictionary of calculated values that you want to compare to the old
+        version. These correspond to the J,H,K,I,V,B magnitudes.
+        
+    Beware: the active isochrone model and the input value isochrone model must be the same for any useful comparison!
+    """
+    ref_libs = json.loads(ref_in.attrs['libraries'])
+    for i in range(2):
+        assert ref_libs[i] == libs[i]
+
+    J1 = np.asarray(ref_in.get('d0/J'))
+    H1 = np.asarray(ref_in.get('d0/H'))
+    K1 = np.asarray(ref_in.get('d0/K'))
+    I1 = np.asarray(ref_in.get('d0/I'))
+    V1 = np.asarray(ref_in.get('d0/V'))
+    B1 = np.asarray(ref_in.get('d0/B'))
+    Ages1 = np.asarray(ref_in.get('t0/Ages'))
+    
+    Jdiff = current['J'] - J1
+    Hdiff = current['H'] - H1
+    Kdiff = current['K'] - K1
+    Idiff = current['I'] - I1
+    Vdiff = current['V'] - V1
+    Bdiff = current['B'] - B1
+    
+    BVdiff = (current['B']-current['V']) - (B1-V1)
+    VIdiff = (current['V']-current['I']) - (V1-I1)
+    VKdiff = (current['V']-current['K']) - (V1-K1)
+    JHdiff = (current['J']-current['H']) - (J1-H1)
+    JKdiff = (current['J']-current['K']) - (J1-K1)
+    HKdiff = (current['H']-current['K']) - (H1-K1)
+
+    print('max B-V difference (input vs. current) = ' + str(np.max(np.abs(BVdiff))))
+    print('max V-I difference (input vs. current) = ' + str(np.max(np.abs(VIdiff))))
+    print('max V-K difference (input vs. current) = ' + str(np.max(np.abs(VKdiff))))
+    print('max J-H difference (input vs. current) = ' + str(np.max(np.abs(JHdiff))))
+    print('max J-K difference (input vs. current) = ' + str(np.max(np.abs(JKdiff))))
+    print('max H-K difference (input vs. current) = ' + str(np.max(np.abs(HKdiff))))
+    
+
+        
+def test_ssp_iso_color(reference_in=None, reference_out=None):
+    """This method will generate a plot of SSP colors as a function of age for
+    different isochrone models, reproducing Conroy & Gunn (2010) Figure 2.
+    We use the MIST, BaSTI, and Padova (the one currently implemented in 
+    python-fsps, not the modified version referenced in Conroy & Gunn (2010))
+    isochrone models. 
+    NOTE: Because changing between different isochrone models requires recompiling FSPS and re-installing python-fsps, we wrote this method to 
+    compute SSP colors for only ONE isochrone model (the active model). The 
+    other two models (the passive models) are provided from a benchmark run on 
+    12/06/2016.
+    The current version of this method assumes that the active model is MIST 
+    and the passive models are BaSTI and Padova. All isochrone models assume 
+    that the spectral library is BaSEL. 
+    NOTE: We leave it to the user to change FSPS libraries to match our active model assumption or to change the active model to suit their needs. Since
+    reference_out outputs only the information from the active model, name 
+    your file appropriately!
+    
+    :param reference_in: (optional)
+        The name (and path) to a file to use as a reference for the quantities
+        being calculated.  The currently calculated quantities will be compared
+        to the ones in this file.
+    :param reference_out: (optional)
+        Name of hdf5 file to output containing the calculated quantities
+    :returns figs:
+        List of matplotlib.pyplot figure objects.
+    """
+    _reset_default_params()
+    
+    # If you want to use a different active model, please change the variable
+    # active_model (in addition to changing the libraries inside FSPS and
+    # python-fsps). Current options are 'MIST', 'Padova', and 'BaSTI'.
+    active_model = 'MIST'
+
+    # Here we will compute SSP J,H,K,I,V,B magnitudes for different ages (and
+    # metallicities)
+    # Conroy & Gunn (2010) only provide metallicities at 7 ages, so we will do
+    # some interpolation - before 1 Gyr, metallicity is constant; after 1 Gyr,
+    # we assume do a linear interpolation
+    met_inter = np.interp(np.linspace(9.1, 10.1, num=21), [9.0, 9.5, 9.8, 10.1], [-0.28, -0.38, -0.68, -1.5])
+    # Full ages and metallicity arrays here
+    logage_act = np.concatenate((np.linspace(7.5, 9.05, num=32), np.linspace(9.1, 10.1, num=21)))
+    met = np.concatenate((np.ones(32)*(-0.28), met_inter))
+    # Magnitude bands
+    bands = ['b','v','cousins_i','2mass_j','2mass_ks','2mass_h']
+    # Initialize some arrays to store the different bands
+    size = len(logage_act)
+    B_act = np.empty(size)
+    V_act = np.empty(size)
+    I_act = np.empty(size)
+    J_act = np.empty(size)
+    K_act = np.empty(size)
+    H_act = np.empty(size)
+    # Now compute the magnitudes
+    print("Computing magnitudes for {} isochrones. This might take a while!".format(active_model))
+    for i in range(size):
+        pop.params['logzsol'] = met[i]
+        mags = pop.get_mags(tage=10**(logage_act[i]-9), bands=bands)
+        B_act[i] = mags[0]
+        V_act[i] = mags[1]
+        I_act[i] = mags[2]
+        J_act[i] = mags[3]
+        K_act[i] = mags[4]
+        H_act[i] = mags[5]
+
+    # Make figures and an output structure
+    # The output structure 
+    fig, axarr = pl.subplots(3,2,figsize=(3.2*2*1.5*0.9, 2.8*3*1.5*0.8))
+    # Colors and linestyles
+    models = {'MIST': {'color':'b', 'linestyle':'-', 'dashes':[8, 4, 2, 4]}, 
+    'Padova': {'color':'k', 'linestyle':'-'}, 
+    'BaSTI':{'color':'r', 'linestyle':'--'}}
+    # Get the passive models data from our benchmark file
+    dir = os.path.dirname(__file__)
+    bench_path = os.path.join(dir, 'benchmark_iso_color.hdf5')
+    f = h5py.File(bench_path, 'r')
+    # Plot all the things!
+    print('Generating plots')
+    for model in models.keys():
+        if model != active_model:
+            # get data for passive models
+            J = np.asarray(f.get(model+'/d0/J'))
+            H = np.asarray(f.get(model+'/d0/H'))
+            K = np.asarray(f.get(model+'/d0/K'))
+            I = np.asarray(f.get(model+'/d0/I'))
+            V = np.asarray(f.get(model+'/d0/V'))
+            B = np.asarray(f.get(model+'/d0/B'))
+            logage = np.asarray(f.get(model+'/t0/Ages'))
+        else:
+            # get data for active model
+            J = J_act
+            H = H_act
+            K = K_act
+            I = I_act
+            V = V_act
+            B = B_act
+            logage = logage_act
+        # plotting
+        axarr[0,0].plot(logage, B-V, lw=2, **models[model])
+        axarr[0,0].set_ylabel(r'${\rm B-V}$', fontsize=13)
+        axarr[0,1].plot(logage, V-I, label='FSPS + '+model, lw=2, \
+            **models[model])
+        axarr[0,1].set_ylabel(r'${\rm V-I}$', fontsize=13)
+        axarr[1,0].plot(logage, V-K, lw=2, **models[model])
+        axarr[1,0].set_ylabel(r'${\rm V-K}$', fontsize=13)
+        axarr[1,1].plot(logage, J-H, lw=2, **models[model])
+        axarr[1,1].set_ylabel(r'${\rm J-H}$', fontsize=13)
+        axarr[2,0].plot(logage, J-K, lw=2, **models[model])
+        axarr[2,0].set_ylabel(r'${\rm J-K}$', fontsize=13)
+        axarr[2,1].plot(logage, H-K, lw=2, **models[model])
+        axarr[2,1].set_ylabel(r'${\rm H-K}$', fontsize=13)
+        
+    #Plot beautification here
+    axarr[0,0].minorticks_on()
+    axarr[0,0].set_xlim(left=7.5,right=10.3)
+    axarr[0,0].set_ylim(top=1.0, bottom=0.0)
+    axarr[0,1].minorticks_on()
+    axarr[0,1].set_xlim(left=7.5,right=10.3)
+    axarr[0,1].set_ylim(top=1.41, bottom=0.3)
+    axarr[1,0].minorticks_on()
+    axarr[1,0].set_xlim(left=7.5,right=10.3)
+    axarr[1,0].set_ylim(top=3.51,bottom=1.01)
+    axarr[1,1].minorticks_on()
+    axarr[1,1].set_xlim(left=7.5,right=10.3)
+    axarr[1,1].set_ylim(top=0.91,bottom=0.35)
+    axarr[2,0].minorticks_on()
+    axarr[2,0].set_xlim(left=7.5,right=10.3)
+    axarr[2,0].set_ylim(top=1.3, bottom=0.39)
+    axarr[2,1].minorticks_on()
+    axarr[2,1].set_xlim(left=7.5,right=10.3)
+    axarr[2,1].set_ylim(top=0.501, bottom=0.0)
+    for i in range(2):
+        axarr[2,i].set_xlabel(r'${\rm log(Age \; [yrs])}$', fontsize=13)
+    leg = axarr[0,1].legend(loc=2, frameon=False)
+    ltext  = leg.get_texts()
+    pl.setp(ltext, fontname='Times New Roman')
+    pl.tight_layout()
+       
+    current = {'B': B_act,
+               'V': V_act,
+               'I': I_act,
+               'J': J_act,
+               'K': K_act,
+               'H': H_act,
+               'Ages': logage_act}
+
+    # Here we compare to an old set of values stored in a benchmark hdf5 file
+    # Please make sure that the log(ages) at which you compute the colors
+    # match our log(ages) array, or else you'll get a lot of errors
+    # We also assume that the input file has data for the magnitudes in the
+    # different bands
+    if reference_in is not None:
+        ref = h5py.File(reference_in, 'r')
+        diffstring = compare_reference_iso(ref, current)
+        ref.close()
+
+    # Here we write out an hdf5 file
+    if reference_out is not None:
+        import datetime
+        timestamp = 'T'.join(str(datetime.datetime.now()).split())
+        ref = h5py.File(reference_out, 'w')
+        # write all the data!
+        ref['d0/J'] = current['J']
+        ref['d0/H'] = current['H']
+        ref['d0/K'] = current['K']
+        ref['d0/I'] = current['I']
+        ref['d0/V'] = current['V']
+        ref['d0/B'] = current['B']
+        ref['d0'].attrs['units'] = 'mags'
+        ref['d0'].attrs['long_name'] = 'Filters J,H,K,I,V,B' 
+        for key in ref['d0'].keys():
+            ref['d0/%s'%key].attrs['units'] = 'mags'
+            ref['d0/%s'%key].attrs['long_name'] = key+'-band magnitude'
+
+        ref['t0/Ages'] = current['Ages']
+        ref['t0'].attrs['units'] = 'Log Age (yrs)'
+        ref['t0'].attrs['long_name'] = 'Log Age of Population in yrs' 
+
+        # now we add useful version things
+        ref.attrs['default'] = 'entry'                        
+        ref.attrs['file_name'] = reference_out
+        ref.attrs['timestamp'] = timestamp
+        ref.attrs['HDF5_Version'] = h5py.version.hdf5_version
+        ref.attrs['h5py_Version'] = h5py.version.version
+        ref.attrs['libraries'] = json.dumps(libs)
+        ref.attrs['default_params'] = json.dumps(default_params)
+        ref.attrs['git_history'] = json.dumps(githashes)
+        ref.close()
+
+    return [fig]

From fba5c5f7280db3d945b3260d3c81425b88aaf9cf Mon Sep 17 00:00:00 2001
From: Irina Ene <irina.ene@berkeley.edu>
Date: Wed, 7 Dec 2016 12:45:23 -0800
Subject: [PATCH 3/4] format some comments

---
 fsps/test_sps_iso.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/fsps/test_sps_iso.py b/fsps/test_sps_iso.py
index fc36dbde..6952a6a6 100644
--- a/fsps/test_sps_iso.py
+++ b/fsps/test_sps_iso.py
@@ -39,7 +39,8 @@ def compare_reference_iso(ref_in, current):
         Dictionary of calculated values that you want to compare to the old
         version. These correspond to the J,H,K,I,V,B magnitudes.
         
-    Beware: the active isochrone model and the input value isochrone model must be the same for any useful comparison!
+    Beware: the active isochrone model and the input value isochrone model
+    must be the same for any useful comparison!
     """
     ref_libs = json.loads(ref_in.attrs['libraries'])
     for i in range(2):
@@ -82,14 +83,16 @@ def test_ssp_iso_color(reference_in=None, reference_out=None):
     We use the MIST, BaSTI, and Padova (the one currently implemented in 
     python-fsps, not the modified version referenced in Conroy & Gunn (2010))
     isochrone models. 
-    NOTE: Because changing between different isochrone models requires recompiling FSPS and re-installing python-fsps, we wrote this method to 
+    NOTE: Because changing between different isochrone models requires
+    recompiling FSPS and re-installing python-fsps, we wrote this method to
     compute SSP colors for only ONE isochrone model (the active model). The 
     other two models (the passive models) are provided from a benchmark run on 
     12/06/2016.
     The current version of this method assumes that the active model is MIST 
     and the passive models are BaSTI and Padova. All isochrone models assume 
     that the spectral library is BaSEL. 
-    NOTE: We leave it to the user to change FSPS libraries to match our active model assumption or to change the active model to suit their needs. Since
+    NOTE: We leave it to the user to change FSPS libraries to match our active
+    model assumption or to change the active model to suit their needs. Since
     reference_out outputs only the information from the active model, name 
     your file appropriately!
     

From 004d96aa27bd5e2130aced9b329062c8250e47e0 Mon Sep 17 00:00:00 2001
From: Irina Ene <irina.ene@berkeley.edu>
Date: Fri, 9 Dec 2016 15:35:09 -0800
Subject: [PATCH 4/4] use arrays to store magnitudes. format some comments.

---
 fsps/test_sps_iso.py | 154 ++++++++++++++++++++-----------------------
 1 file changed, 73 insertions(+), 81 deletions(-)

diff --git a/fsps/test_sps_iso.py b/fsps/test_sps_iso.py
index 6952a6a6..7926aa0e 100644
--- a/fsps/test_sps_iso.py
+++ b/fsps/test_sps_iso.py
@@ -37,37 +37,45 @@ def compare_reference_iso(ref_in, current):
         h5py.File object with keys matching those in the given dictionary of values.
     :param current:
         Dictionary of calculated values that you want to compare to the old
-        version. These correspond to the J,H,K,I,V,B magnitudes.
+        version. These correspond to the J,H,K,I,V,B magnitudes and the ages
+        at which the magnitudes are computed.
         
-    Beware: the active isochrone model and the input value isochrone model
-    must be the same for any useful comparison!
+    Beware: the active isochrone model and the input value isochrone model must be the same for any useful comparison!
     """
     ref_libs = json.loads(ref_in.attrs['libraries'])
     for i in range(2):
         assert ref_libs[i] == libs[i]
 
-    J1 = np.asarray(ref_in.get('d0/J'))
-    H1 = np.asarray(ref_in.get('d0/H'))
-    K1 = np.asarray(ref_in.get('d0/K'))
-    I1 = np.asarray(ref_in.get('d0/I'))
-    V1 = np.asarray(ref_in.get('d0/V'))
-    B1 = np.asarray(ref_in.get('d0/B'))
-    Ages1 = np.asarray(ref_in.get('t0/Ages'))
-    
-    Jdiff = current['J'] - J1
-    Hdiff = current['H'] - H1
-    Kdiff = current['K'] - K1
-    Idiff = current['I'] - I1
-    Vdiff = current['V'] - V1
-    Bdiff = current['B'] - B1
+    # create empty arrays to store data
+    data_in = np.zeros((len(ref_in['d0'].keys())+1, \
+        len(np.asarray(ref_in.get('t0/Ages')))))
+    data_c = np.zeros((len(current.keys())+1, \
+        len(current['Ages'])))
+
+    # populate with ages
+    data_in[0,:] = np.asarray(ref_in.get('t0/Ages'))
+    data_c[0,:] = current['Ages']
     
-    BVdiff = (current['B']-current['V']) - (B1-V1)
-    VIdiff = (current['V']-current['I']) - (V1-I1)
-    VKdiff = (current['V']-current['K']) - (V1-K1)
-    JHdiff = (current['J']-current['H']) - (J1-H1)
-    JKdiff = (current['J']-current['K']) - (J1-K1)
-    HKdiff = (current['H']-current['K']) - (H1-K1)
+    # sort keys
+    sorted_keys = current.keys()
+    sorted_keys.remove('Ages')
+    sorted_keys.sort()
+    # populate with magnitudes
+    for i, key in enumerate(sorted_keys):
+        data_in[i+1,:] = np.asarray(ref_in.get('d0/%s'%key))
+        data_c[i+1,:] = current[key]
+
+    # compute color difference between current and input
+    BVdiff = (data_c[1,:]-data_c[6,:]) - (data_in[1,:]-data_in[6,:])
+    VIdiff = (data_c[6,:]-data_c[3,:]) - (data_in[6,:]-data_in[3,:])
+    VKdiff = (data_c[6,:]-data_c[5,:]) - (data_in[6,:]-data_in[5,:])
+    JHdiff = (data_c[4,:]-data_c[2,:]) - (data_in[4,:]-data_in[2,:])
+    JKdiff = (data_c[4,:]-data_c[5,:]) - (data_in[4,:]-data_in[5,:])
+    HKdiff = (data_c[2,:]-data_c[5,:]) - (data_in[2,:]-data_in[5,:])
 
+    #colordiffset = np.asarray(BVdiff, VIdiff, VKdiff, JHdiff, JKdiff, HKdiff)
+    
+    
     print('max B-V difference (input vs. current) = ' + str(np.max(np.abs(BVdiff))))
     print('max V-I difference (input vs. current) = ' + str(np.max(np.abs(VIdiff))))
     print('max V-K difference (input vs. current) = ' + str(np.max(np.abs(VKdiff))))
@@ -75,7 +83,7 @@ def compare_reference_iso(ref_in, current):
     print('max J-K difference (input vs. current) = ' + str(np.max(np.abs(JKdiff))))
     print('max H-K difference (input vs. current) = ' + str(np.max(np.abs(HKdiff))))
     
-
+    #return data_in, data_c, diffset, colordiffset
         
 def test_ssp_iso_color(reference_in=None, reference_out=None):
     """This method will generate a plot of SSP colors as a function of age for
@@ -83,16 +91,14 @@ def test_ssp_iso_color(reference_in=None, reference_out=None):
     We use the MIST, BaSTI, and Padova (the one currently implemented in 
     python-fsps, not the modified version referenced in Conroy & Gunn (2010))
     isochrone models. 
-    NOTE: Because changing between different isochrone models requires
-    recompiling FSPS and re-installing python-fsps, we wrote this method to
+    NOTE: Because changing between different isochrone models requires recompiling FSPS and re-installing python-fsps, we wrote this method to 
     compute SSP colors for only ONE isochrone model (the active model). The 
     other two models (the passive models) are provided from a benchmark run on 
     12/06/2016.
     The current version of this method assumes that the active model is MIST 
     and the passive models are BaSTI and Padova. All isochrone models assume 
     that the spectral library is BaSEL. 
-    NOTE: We leave it to the user to change FSPS libraries to match our active
-    model assumption or to change the active model to suit their needs. Since
+    NOTE: We leave it to the user to change FSPS libraries to match our active model assumption or to change the active model to suit their needs. Since
     reference_out outputs only the information from the active model, name 
     your file appropriately!
     
@@ -117,39 +123,30 @@ def test_ssp_iso_color(reference_in=None, reference_out=None):
     # Conroy & Gunn (2010) only provide metallicities at 7 ages, so we will do
     # some interpolation - before 1 Gyr, metallicity is constant; after 1 Gyr,
     # we assume do a linear interpolation
-    met_inter = np.interp(np.linspace(9.1, 10.1, num=21), [9.0, 9.5, 9.8, 10.1], [-0.28, -0.38, -0.68, -1.5])
+    met_inter = np.interp(np.linspace(9.1, 10.1, num=21), \
+        [9.0, 9.5, 9.8, 10.1], [-0.28, -0.38, -0.68, -1.5])
     # Full ages and metallicity arrays here
-    logage_act = np.concatenate((np.linspace(7.5, 9.05, num=32), np.linspace(9.1, 10.1, num=21)))
+    logage_act = np.concatenate((np.linspace(7.5, 9.05, num=32), \
+        np.linspace(9.1, 10.1, num=21)))
     met = np.concatenate((np.ones(32)*(-0.28), met_inter))
     # Magnitude bands
     bands = ['b','v','cousins_i','2mass_j','2mass_ks','2mass_h']
     # Initialize some arrays to store the different bands
-    size = len(logage_act)
-    B_act = np.empty(size)
-    V_act = np.empty(size)
-    I_act = np.empty(size)
-    J_act = np.empty(size)
-    K_act = np.empty(size)
-    H_act = np.empty(size)
+    t_size = len(logage_act)
+    colors_act = np.empty([len(bands), t_size])
     # Now compute the magnitudes
     print("Computing magnitudes for {} isochrones. This might take a while!".format(active_model))
-    for i in range(size):
+    for i in range(t_size):
         pop.params['logzsol'] = met[i]
-        mags = pop.get_mags(tage=10**(logage_act[i]-9), bands=bands)
-        B_act[i] = mags[0]
-        V_act[i] = mags[1]
-        I_act[i] = mags[2]
-        J_act[i] = mags[3]
-        K_act[i] = mags[4]
-        H_act[i] = mags[5]
+        colors_act[:,i] = pop.get_mags(tage=10**(logage_act[i]-9), bands=bands)
 
     # Make figures and an output structure
     # The output structure 
     fig, axarr = pl.subplots(3,2,figsize=(3.2*2*1.5*0.9, 2.8*3*1.5*0.8))
     # Colors and linestyles
     models = {'MIST': {'color':'b', 'linestyle':'-', 'dashes':[8, 4, 2, 4]}, 
-    'Padova': {'color':'k', 'linestyle':'-'}, 
-    'BaSTI':{'color':'r', 'linestyle':'--'}}
+              'Padova': {'color':'k', 'linestyle':'-'}, 
+              'BaSTI':{'color':'r', 'linestyle':'--'}}
     # Get the passive models data from our benchmark file
     dir = os.path.dirname(__file__)
     bench_path = os.path.join(dir, 'benchmark_iso_color.hdf5')
@@ -157,37 +154,35 @@ def test_ssp_iso_color(reference_in=None, reference_out=None):
     # Plot all the things!
     print('Generating plots')
     for model in models.keys():
+        temp_plot = np.empty([len(bands)+1, t_size])
         if model != active_model:
             # get data for passive models
-            J = np.asarray(f.get(model+'/d0/J'))
-            H = np.asarray(f.get(model+'/d0/H'))
-            K = np.asarray(f.get(model+'/d0/K'))
-            I = np.asarray(f.get(model+'/d0/I'))
-            V = np.asarray(f.get(model+'/d0/V'))
-            B = np.asarray(f.get(model+'/d0/B'))
-            logage = np.asarray(f.get(model+'/t0/Ages'))
+            temp_plot[0,:] = np.asarray(f.get(model+'/t0/Ages'))
+            for i, key in enumerate(['B', 'V', 'I', 'J', 'K', 'H']):
+                temp_plot[i+1,:] = np.asarray(f.get(model+'/d0/%s'%key))
         else:
             # get data for active model
-            J = J_act
-            H = H_act
-            K = K_act
-            I = I_act
-            V = V_act
-            B = B_act
-            logage = logage_act
+            temp_plot[0,:] = logage_act
+            for i in range(len(bands)):
+                temp_plot[i+1,:] = colors_act[i,:]
         # plotting
-        axarr[0,0].plot(logage, B-V, lw=2, **models[model])
-        axarr[0,0].set_ylabel(r'${\rm B-V}$', fontsize=13)
-        axarr[0,1].plot(logage, V-I, label='FSPS + '+model, lw=2, \
+        axarr[0,0].plot(temp_plot[0,:], temp_plot[1,:]-temp_plot[2,:], lw=2, \
             **models[model])
+        axarr[0,0].set_ylabel(r'${\rm B-V}$', fontsize=13)
+        axarr[0,1].plot(temp_plot[0,:], temp_plot[2,:]-temp_plot[3,:], \
+            label='FSPS + {}'.format(model), lw=2, **models[model])
         axarr[0,1].set_ylabel(r'${\rm V-I}$', fontsize=13)
-        axarr[1,0].plot(logage, V-K, lw=2, **models[model])
+        axarr[1,0].plot(temp_plot[0,:], temp_plot[2,:]-temp_plot[5,:], lw=2, \
+            **models[model])
         axarr[1,0].set_ylabel(r'${\rm V-K}$', fontsize=13)
-        axarr[1,1].plot(logage, J-H, lw=2, **models[model])
+        axarr[1,1].plot(temp_plot[0,:], temp_plot[4,:]-temp_plot[6,:], lw=2, \
+            **models[model])
         axarr[1,1].set_ylabel(r'${\rm J-H}$', fontsize=13)
-        axarr[2,0].plot(logage, J-K, lw=2, **models[model])
+        axarr[2,0].plot(temp_plot[0,:], temp_plot[4,:]-temp_plot[5,:], lw=2, \
+            **models[model])
         axarr[2,0].set_ylabel(r'${\rm J-K}$', fontsize=13)
-        axarr[2,1].plot(logage, H-K, lw=2, **models[model])
+        axarr[2,1].plot(temp_plot[0,:], temp_plot[6,:]-temp_plot[5,:], lw=2, \
+            **models[model])
         axarr[2,1].set_ylabel(r'${\rm H-K}$', fontsize=13)
         
     #Plot beautification here
@@ -216,12 +211,12 @@ def test_ssp_iso_color(reference_in=None, reference_out=None):
     pl.setp(ltext, fontname='Times New Roman')
     pl.tight_layout()
        
-    current = {'B': B_act,
-               'V': V_act,
-               'I': I_act,
-               'J': J_act,
-               'K': K_act,
-               'H': H_act,
+    current = {'B': colors_act[0,:],
+               'V': colors_act[1,:],
+               'I': colors_act[2,:],
+               'J': colors_act[3,:],
+               'K': colors_act[4,:],
+               'H': colors_act[5,:],
                'Ages': logage_act}
 
     # Here we compare to an old set of values stored in a benchmark hdf5 file
@@ -240,14 +235,11 @@ def test_ssp_iso_color(reference_in=None, reference_out=None):
         timestamp = 'T'.join(str(datetime.datetime.now()).split())
         ref = h5py.File(reference_out, 'w')
         # write all the data!
-        ref['d0/J'] = current['J']
-        ref['d0/H'] = current['H']
-        ref['d0/K'] = current['K']
-        ref['d0/I'] = current['I']
-        ref['d0/V'] = current['V']
-        ref['d0/B'] = current['B']
+        for key in ['B', 'V', 'I', 'J', 'K', 'H']:
+            ref['d0/'+key] = current[key]
+            
         ref['d0'].attrs['units'] = 'mags'
-        ref['d0'].attrs['long_name'] = 'Filters J,H,K,I,V,B' 
+        ref['d0'].attrs['long_name'] = 'Filters B,V,I,J,K,H' 
         for key in ref['d0'].keys():
             ref['d0/%s'%key].attrs['units'] = 'mags'
             ref['d0/%s'%key].attrs['long_name'] = key+'-band magnitude'