Created
April 19, 2017 15:41
Revisions
-
ycopin created this gist
Apr 19, 2017 .There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,226 @@ # lbda: 4327A # radius[spx] flux 0.716608613055 12.3552 0.745272338424 11.5852 0.745343175736 11.7769 0.834997936881 10.4475 1.55358944456 4.54829 1.57210666273 4.74318 1.61169010512 4.2645 1.68493686848 3.5319 1.74032685753 3.35422 1.75684694171 3.12297 1.79226673209 3.36554 1.85843829222 3.02449 2.22949482886 1.52675 2.22956586929 1.63293 2.26810158109 1.49356 2.3864965395 1.25747 2.50436608068 0.936273 2.53564800666 1.00553 2.57673311389 1.15053 2.64212091831 0.831548 2.84900419776 0.594853 2.8765223231 0.583048 2.91271272972 0.675153 2.95450069299 0.464054 2.97073230766 0.661472 2.98549666889 0.826283 3.03348947381 0.520949 3.15229100498 0.431523 3.15583097812 0.429738 3.18473553815 0.40955 3.22985083738 0.390899 3.34164711507 0.400866 3.48282339057 0.264776 3.51958836633 0.289079 3.56111809486 0.369776 3.62251206252 0.301754 3.71499467699 0.209787 3.71506573374 0.280393 3.80624506539 0.181903 3.81961937739 0.245014 3.85427253245 0.315292 3.90638349286 0.21185 3.93801478407 0.216523 3.97741102894 0.182418 4.00963045684 0.190295 4.02225203426 0.201376 4.04604214996 0.215724 4.10019441344 0.223117 4.19075895105 0.159204 4.23430276792 0.167069 4.28187742792 0.166068 4.38779345193 0.172258 4.41132978275 0.119673 4.43217301666 0.160119 4.55231486508 0.163044 4.55489274787 0.138129 4.61133216536 0.156659 4.61574725029 0.130903 4.63553478918 0.115133 4.67432513177 0.132926 4.75308449483 0.135578 4.77468775407 0.145954 4.86763479788 0.12417 4.94166816714 0.119244 5.11250019988 0.100622 5.14730787421 0.101647 5.18381181554 0.121259 5.20067742058 0.0851084 5.20074848183 0.110162 5.23572576045 0.120342 5.24276662586 0.113552 5.27072464047 0.0976773 5.3226859422 0.101061 5.37114048618 0.102004 5.37139631847 0.0972523 5.47300132818 0.102482 5.4895360163 0.0946428 5.494603756 0.0905185 5.54667280234 0.0987886 5.60412125453 0.104601 5.60565947748 0.0864304 5.63826802659 0.0823288 5.72187640056 0.087733 5.76272665222 0.0884338 5.87321571525 0.087444 5.88741514768 0.090112 5.89834272182 0.0884296 6.08820976309 0.0842796 6.10320793027 0.0701937 6.10975008792 0.0780519 6.2109695558 0.071304 6.25063231368 0.0781598 6.28707456055 0.0712531 6.29016204293 0.076421 6.32362546528 0.0949923 6.36647668914 0.0752858 6.37287235961 0.0686733 6.37408350084 0.0894448 6.40590213299 0.0689079 6.43020509584 0.0848073 6.47952226162 0.071904 6.54275077635 0.0773923 6.57796639746 0.0802376 6.59908825003 0.0705312 6.64504325816 0.0650343 6.67331319186 0.0732772 6.6844667654 0.0713967 6.68631708223 0.0718472 6.68649220901 0.069661 6.81477152573 0.0720635 6.8485932839 0.0717355 6.92172979656 0.07091 7.0012341749 0.0654785 7.01223861405 0.0660574 7.04105826117 0.0682567 7.04980125208 0.0667779 7.07574269414 0.0687417 7.27140842143 0.0642235 7.27310203124 0.0585062 7.37238482467 0.0664186 7.38586791292 0.0680497 7.39040165176 0.0707092 7.42796882965 0.0668018 7.48741291865 0.0565187 7.50171787322 0.0633702 7.53986712527 0.0627339 7.56589687247 0.0628192 7.59537730648 0.0637952 7.59811853646 0.0637862 7.65665216215 0.0593574 7.71392687641 0.0592076 7.75354246372 0.0620853 7.81951423141 0.0661176 7.83954189482 0.0555423 7.89228038364 0.0632143 7.93813168346 0.0649871 7.94932363649 0.0656511 8.05903117593 0.0696435 8.12372685071 0.0631257 8.17226365677 0.064036 8.19845933265 0.0591071 8.29727509571 0.0531699 8.36770178377 0.0604193 8.50830765398 0.0608632 8.52398615886 0.055592 8.53115203558 0.0583338 8.53765804024 0.0594891 8.56954033569 0.0563404 8.59252856405 0.0573073 8.59258097012 0.0594579 8.60953474605 0.0597631 8.63519245452 0.0559289 8.71203502457 0.0545608 8.72315814804 0.0583596 8.79695263499 0.0555135 8.80156929375 0.0556808 8.85036671934 0.0563155 8.88011142202 0.057412 8.92400234686 0.0583789 8.98609785093 0.0564457 9.04624865759 0.0603954 9.10006349 0.0585736 9.14090205762 0.0555651 9.26072865572 0.0593005 9.29909399773 0.0577292 9.33805436566 0.0581372 9.34914012608 0.0560803 9.49498230705 0.0562901 9.53591177379 0.0545264 9.58499289816 0.0566096 9.59027308371 0.0529429 9.60940640301 0.0549165 9.61814821684 0.0503754 9.65805003171 0.063911 9.67253636646 0.0561432 9.71155025843 0.0575919 9.73661702538 0.0582254 9.76118321868 0.0579715 9.77895436687 0.0560509 9.83355597488 0.0562717 9.85243406899 0.0592493 9.9017908702 0.0534869 9.97487957323 0.0558348 10.0555383063 0.0562221 10.0909219472 0.0512105 10.0944264535 0.0565847 10.163130117 0.0565575 10.1736796356 0.0571313 10.1852922705 0.0542905 10.3305975653 0.0535239 10.3641298066 0.0529833 10.4203145061 0.0548901 10.5134658945 0.0561336 10.5252947021 0.0495351 10.528141146 0.055008 10.534496756 0.054671 10.5884431583 0.0561611 10.6042916555 0.0506579 10.6265506055 0.0567964 10.764270591 0.0507111 10.7954280765 0.0549711 10.8514769538 0.0568192 10.9442549519 0.0561116 11.0528078648 0.0602084 11.0561657198 0.0491225 11.0580784016 0.0556615 11.1016403932 0.0562249 11.1438453632 0.0571678 11.1784741553 0.0519887 11.4464172469 0.0586112 11.4552593037 0.0509309 11.5137508099 0.0533362 11.6788378502 0.0546892 11.8122319215 0.051182 11.8155777402 0.0591309 11.9508165928 0.0568809 12.011004085 0.0511922 12.0208727733 0.052377 12.5248695186 0.055433 12.5313987639 0.0574275 12.6115499831 0.0611926 12.6296464902 0.0517462 13.2746333263 0.0621648 13.2948005553 0.0534115 This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,226 @@ # lbda: 7603A # radius[spx] flux 0.418593273815 18.1837 0.5349266477 16.3112 0.919587991512 9.29802 0.971844008792 9.15756 1.22670232071 5.66755 1.27585746046 5.48179 1.28633172146 7.10789 1.40976163455 5.59343 1.52877680815 3.99235 1.62286001831 3.41529 1.72398957832 2.70682 1.82788979249 2.72601 1.87511701485 2.20257 1.89810827837 2.46238 2.17357618332 1.56991 2.19872163731 1.40708 2.22923260348 1.38511 2.24046150617 1.12481 2.29762145977 0.996929 2.29780795798 1.19263 2.3278810413 1.03583 2.43171208388 0.821884 2.4555725823 0.93069 2.50791364179 0.71838 2.57777134088 0.674734 2.58745711476 0.741269 2.84844288343 0.460966 2.86120351681 0.623048 2.8742215112 0.704532 2.9235907623 0.46231 3.05536249719 0.34627 3.06342499664 0.3748 3.11483177589 0.404891 3.15058219771 0.453676 3.16397626403 0.330685 3.17667387593 0.318302 3.1876411399 0.272309 3.18809584127 0.37018 3.21109696855 0.404542 3.26770385838 0.2804 3.28376732383 0.237104 3.39690804064 0.197484 3.40692260011 0.254243 3.47112194369 0.296218 3.57742101492 0.184147 3.57885785942 0.225812 3.63063936958 0.214682 3.661311191 0.168465 3.73780741908 0.169168 3.82630935387 0.162878 3.85446589165 0.252423 3.86224189251 0.136527 3.86257473124 0.272684 3.92129908485 0.171978 3.95412027557 0.128086 4.04398798761 0.129368 4.04780017091 0.185285 4.07092310874 0.185173 4.07840690215 0.112517 4.11294243377 0.111917 4.1525157644 0.101422 4.18252714849 0.162912 4.18385292786 0.104878 4.20152466164 0.178031 4.24539969208 0.0883317 4.24585344043 0.106553 4.35024279457 0.111205 4.36722000066 0.0891081 4.38913660053 0.092006 4.40482017513 0.131164 4.41619975722 0.109993 4.45148306955 0.115698 4.45571223124 0.0770115 4.4941655258 0.0776948 4.5049795378 0.0874144 4.62746520586 0.0598339 4.68466077254 0.0774905 4.76890209863 0.0870922 4.84519522369 0.0633349 4.85049781519 0.124598 4.85569951845 0.113245 4.8776410512 0.0595956 4.92003820434 0.0642767 4.92049891027 0.0606985 4.94604446013 0.0569694 4.96481049366 0.0716581 4.96951756618 0.0603814 4.97444042875 0.0567405 5.00671503634 0.095196 5.02181856802 0.0877579 5.0293142344 0.0584706 5.04764817346 0.0530278 5.10864052491 0.0463022 5.15342952743 0.0523819 5.17720483634 0.0427152 5.17910335895 0.0883288 5.19014788889 0.0477646 5.19562077363 0.0913007 5.26049131125 0.0565143 5.28127811199 0.0398706 5.29284407504 0.0404179 5.31012401111 0.0636266 5.31432998036 0.0629292 5.33384109296 0.0632562 5.36248020998 0.0714763 5.43141889242 0.034259 5.52359972433 0.0344777 5.57434660886 0.0406432 5.59029961952 0.0419796 5.62747672452 0.0338062 5.66641579305 0.0514407 5.69998204593 0.0343423 5.72363349568 0.0342854 5.73647291599 0.0393544 5.73742047701 0.0396858 5.76813791775 0.0476922 5.76917145059 0.0306585 5.8008678815 0.0340075 5.84501281891 0.0340159 5.84751790108 0.0321971 5.84788430121 0.0645667 5.85116672745 0.0579969 5.8608158118 0.0367026 5.90782653411 0.0246767 5.91058130148 0.0369088 5.97910567615 0.0552103 5.98024687102 0.0267481 5.98758312595 0.0314273 5.9887317942 0.0479029 6.03266233444 0.0287829 6.08692837443 0.0308964 6.11728928452 0.0245978 6.12542085599 0.0268558 6.17678636974 0.0550516 6.19161860618 0.0504805 6.19672118636 0.0247216 6.23634618436 0.0381883 6.25172137027 0.0363088 6.2520150473 0.0280512 6.2591966974 0.0238032 6.2632998974 0.023682 6.28964179494 0.0360326 6.29947017055 0.0297034 6.33324144729 0.0429346 6.36808291794 0.0200404 6.41373226008 0.0201326 6.48391923936 0.0206196 6.48895908482 0.0215115 6.50419230415 0.0190411 6.52356138304 0.0229855 6.53286166452 0.0199665 6.59348986905 0.0288672 6.60177864589 0.0196494 6.60489802601 0.0279993 6.62521598429 0.0294355 6.6278715614 0.0281915 6.67770424287 0.0174782 6.68458791426 0.0239929 6.7211343749 0.0259512 6.86514476453 0.015365 6.86618289053 0.0177741 6.8875793812 0.0197852 6.90547528991 0.0195681 6.91945866139 0.017058 6.95834844892 0.0172877 6.9590678424 0.0180712 6.97799393422 0.0210625 7.05844323299 0.0148512 7.06649381696 0.0149227 7.06736848445 0.0188029 7.09177754472 0.022784 7.17511450167 0.0321384 7.1887277529 0.031752 7.1907315832 0.0145086 7.21151333727 0.0120302 7.22452655067 0.01683 7.27166313007 0.0226572 7.3022213752 0.0140301 7.31188761244 0.0245256 7.32554291419 0.0143192 7.41441143823 0.0121539 7.47410572581 0.0132615 7.52474948616 0.0119247 7.539227372 0.0166344 7.56823682819 0.0176403 7.60664655512 0.0114279 7.63436217423 0.0187077 7.67781336256 0.0128527 7.69365896062 0.0111624 7.77417412219 0.0105209 7.80439266302 0.013102 7.81605793843 0.0109135 7.83709865142 0.0102302 7.86171440812 0.0100137 7.94409994108 0.00996457 8.16110686869 0.00901483 8.17385137906 0.021416 8.18654215444 0.0209207 8.20760631927 0.00945191 8.22823505858 0.0101361 8.25800069301 0.0142736 8.26822306553 0.00766784 8.29562858885 0.0145939 8.37576925958 0.00988022 8.3934165594 0.00899557 8.43609267865 0.00797652 8.49739901854 0.0104597 8.53430172484 0.00977896 8.62321880277 0.00827871 8.70236169544 0.00597591 8.74704676139 0.00694602 8.78546991541 0.00552142 8.85770209073 0.00976665 9.0221104261 0.0074013 9.0490272202 0.0057893 9.14912536131 0.0073833 9.15169538104 0.00507639 9.28700686923 0.00642731 9.46726166027 0.00552872 9.58712314609 0.00491504 9.72848613052 0.00694392 9.94804243126 0.00547246 10.0827916396 0.00495822 10.4827672035 0.00603948 This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal file line number Diff line number Diff line change @@ -0,0 +1,545 @@ #/usr/bin/env python # -*- coding: utf-8 -*- """ Simulation of seeing in large telescopes (D>>r0), tested against SNfactory ad-hoc PSF models. References: * Tokovinin 2002PASP..114.1156T * Fried 1966JOSA...56.1372F * Racine 1996PASP..108..699R Author: Yannick Copin <y.copin@ipnl.in2p3.fr> """ from __future__ import division import numpy as N from scipy import special as S from scipy import fftpack as F from scipy import interpolate as I from scipy.optimize import fsolve from scipy.integrate import quad import matplotlib.pyplot as P import hankel as H import fractions RAD2ARCS = 206264.80624709636 # From radians to arcsecs # Color names BLUE = 'C0' GREEN = 'C2' RED = 'C3' def otf_Kolmogorov(r, r0=0.15, expo=5/3): """ Kolmogorov optical transfer function. r [m/rad], r0 is the Fried parameter [m]. Note: lower exponent (e.g. 3/2=1.5 instead of 5/3=1.67) means more wings. 1.55 corresponds roughly to long exposure ad-hoc SNfactory PSFs, short SNfactory PSFs are closer to 1.5. """ # (Half) Phase structure function hdphi = (r/r0)**expo if expo == (5/3): # 2*( 24/5*gamma(6/5) )**(5/6) = 6.8838771822938112... # Including the 0.5 of otf: 3.4419385911469056... hdphi *= 3.4419385911469056 else: # Possible generalisation, including the 0.5 of otf hdphi *= (8/expo*S.gamma(2/expo))**(expo/2) # otf = exp(-0.5*phase), the 0.5 factor has already been included return N.exp(-hdphi) def otf_VonKarman(r, r0=0.15, L0=20., desaturate=True): """ Von Karman optical transfer function. Note that Tokovinin's formula has a typo -- 2**(11/6) should be 2**(5/6) --, see Conan (2000). """ #raise NotImplementedError("Von Karman phase does not work properly...") print "WARNING: Von Karman phase does not work properly..." # (Half) Phase structure function def hphase(r): x = 2*N.pi*r/L0 # gamma(11/6)/(2**(5/6)*pi**(8/3)) * ( 24/5*gamma(6/5) )**(5/6) = # 0.08583068106228546... # gamma(5/6)/2**(1/6) = 1.0056349179985902... return 0.08583068106228546 * (r0/L0)**(-5/3) * \ ( 1.0056349179985902 - x**(5/6) * S.kv(5/6, x) ) hdphi = hphase(r) hdphi[r == 0] = 0. # Extension to 0 (NaN otherwise) otf = N.exp(-hdphi) # Optical transfert function if desaturate: # See Tokovinin (2002) otf -= N.exp(-hphase(N.max(r)*2)) # Remove constant level otf /= otf.max() # Renormalization return otf def psf2D_FFT(otfarr, normed=True): """ FFT-computed 2D-PSF from input OTF 2D-array (e.g. from otf_Kolmogorov). """ assert otfarr.ndim == 2, "Input OTF array is not 2D" psf = F.fftshift(N.absolute(F.fft2(otfarr))) # PSF=FT(OTF), otfarr.shape if normed: psf /= psf.max() return psf def friedParamater(lbda, r0ref=0.10, lref=5e-7, expo=5/3): """Implement r0 chromaticity.""" return r0ref*(lbda/lref)**(2/expo) def seeing_fwhm(lbda, r0ref=0.10, lref=5e-7, expo=5/3, verbose=False): """ Seeing FWHM [arcsec] at wavelength lbda [m] for Fried parameter r0ref [m] at reference wavelength lref [m]. Note: the 'traditional' seeing is dependant on exponent 5/3 used in otf_Kolmogorov through constant 0.976 (actually 0.97586...). """ r0 = friedParamater(lbda, r0ref, lref) s = 0.975863432579*lbda/r0 * RAD2ARCS # [arcsec] Valid for expo=5/3 if expo != (5/3): s = 2*fsolve(lambda r: psf_Kolmogorov_Hankel(r, lbda, r0=r0, expo=expo) - 0.5, s/2) if verbose: print "r0 = %.0f cm @%.0fÅ = %.0f cm @%.0fÅ: " \ "seeing=%.2f\" (FWHM) (expo=%.2f)" % \ (r0ref*100, lref*1e10, r0*100, lbda*1e10, s, expo) return s def psf_ES(r, alpha, coeffs, decompose=False): """ Compute (max-normalized) extract_star PSF given free parameter alpha and set of correlation coefficients coeffs for radius array r [arcsec]. If decompose=True, return individual components of PSF (gaussian + moffat). """ spx = 0.43 # Spaxel size [arcsec] name, s1, s0, b1, b0, e1, e0 = coeffs sigma = s0 + s1*alpha beta = b0 + b1*alpha eta = e0 + e1*alpha print "PSF %s: sigma,beta,eta =" % name, sigma, beta, eta r2 = (r/spx)**2 # (Radius [spx])**2 gaussian = N.exp(-r2/2/sigma**2) moffat = (1 + r2/alpha**2)**(-beta) s = eta*gaussian + moffat if decompose: return eta*gaussian/s.max(), moffat/s.max() else: return s/s.max() def psf_ES_long(r, alpha, decompose=False): return psf_ES(r, alpha, ['long', 0.215, 0.545, 0.345, 1.685, 0.0, 1.04], decompose=decompose) def psf_ES_short(r, alpha, decompose=False): return psf_ES(r, alpha, ['short', 0.2, 0.56, 0.415, 1.395, 0.16, 0.6], decompose=decompose) def psf_Kolmogorov_FFT(r, lbda, scale=None, r0=0.15, expo=5/3): """ Compute (max-normalized) Kolmogorov (FFT-based) PSF given free parameters r0 and expo for radius array r [arcsec] at wavelength lbda [m]. """ assert scale is not None, "Pixel scale [arcsec/px] should be set" amax = N.max(r) # Image half width [arcsec] if amax < 20: # print "WARNING: extending PSF modeling from %.2f'' to 20''" % amax amax = 20 # 2D-array construction xmax = int(amax / scale) # Image half width [px] amax = xmax * scale # Image half width [arcsec] # print "Image half width: %d px = %.2f\"" % (xmax,amax) x = N.arange(-xmax, xmax+0.1, dtype='d') # (n,) [px] rpx = N.hypot(*N.meshgrid(x, x)) # Radius 2D-array (n,n) [px] px2f = RAD2ARCS/scale/len(x) # [px/rad] f = rpx * px2f # Spatial freq 2D-array (n,n) [1/rad] # FFT-based PSF computation psf = psf2D_FFT(otf_Kolmogorov(lbda*f, r0=r0, expo=expo)) # (n,n) a = rpx * scale # Radius (n,n) [arcsec] j = N.argsort(a.ravel()) # Such that a.ravel()[j] is sorted i = j[:N.searchsorted(a.ravel()[j], amax)] x = a.ravel()[i] y = N.log(psf.ravel()[i]) # Spline interpolation n, = N.nonzero(N.diff(x)) # Discard duplicated radii (and last point) spl = I.UnivariateSpline(x[n], y[n], k=3, s=0) # Estimate PSF at given radii return N.exp(spl(r)) def psf_Kolmogorov_Hankel(r, lbda, r0=0.15, expo=5/3, normed=True): """ Compute (max-normalized) Kolmogorov (Hankel-based) PSF given free parameters r0 and expo for radius array r [arcsec] at wavelength lbda [m]. """ # One should compute psf(r) = H.hankel(lambda # r:otf_Kolmogorov(lbda*r/(2*N.pi)), r/RAD2ARCS), but for some numerical # reasons, one computes lbda**2/(2*pi) * psf(r) = H.hankel(lambda # r:otf_Kolmogorov(r, r0=r0, expo=expo), 2*pi/lambda*r/RAD2ARCS) tmp = (6.2831853071795862/RAD2ARCS/lbda)*N.atleast_1d(r) psf = H.hankel(lambda r: otf_Kolmogorov(r, r0=r0, expo=expo), tmp) # psf(0) = 2*pi/expo*fO**2*Gamma(2/expo) # where f0 = r0/lbda * beta**(-1/expo) # and beta = (8/expo*Gamma(2/expo))**(expo/2) = 3.44... for expo=5/3 # Therefore f0**2 = (r0/lbda)**2 / ( 8/expo*Gamma(2/expo) ) # and psf(0) = pi/4 * (r0/lbda)**2. Given we actually compute # lbda**2/(2*pi)*psf(r), this gives psf(0) = r0**2/8 pmax = r0**2/8 psf[tmp == 0] = pmax # Remove NaNs at r=0 if normed: psf /= pmax # Max normalisation return psf.reshape(N.shape(r)) # PSF should have the same shape as r def integ_radial_laguerre(f, accuracy=1e-3, imin=1, imax=40): """ Compute s = 2*pi int_0^infinity g(r) dr with g(r) = r*f(r), using Gauss-Laguerre quadrature: int_0^infinity exp(-r)*exp(r)*g(r) dr = sum_i w_i exp(r_i) g(r_i) where r_i and w_i are the roots and weights of ith-order Laguerre polynomial. """ def g(r): return r*f(r) sprev = 0 errprev = N.inf for i in xrange(imin, imax): w = S.laguerre(i+1).weights # Laguerre polynomial roots & weights s = N.dot(w[:, 2], g(w[:, 0])) # sum_i w_i exp(x_i) g(x_i) err = abs(1-sprev/s) print "Order %d (%f-%f): s=%f err=%f" % \ (i+1, min(w[:, 0]), max(w[:, 0]), 2*N.pi*s, err) if err < accuracy: print "Requested accuracy is reached" break elif err > errprev: print "Relative error is increasing" # break sprev = s errprev = err s *= 6.2831853071795862 # 2*pi return s, s*err # Integral and error estimate def integ_radial_quad(f, epsabs=1e-6, epsrel=1e-6): """Compute s = 2*pi int_0^infinity g(r) dr with g(r) = r*f(r).""" s, ds = quad(lambda r: r*f(r), 0, N.inf, epsabs=epsabs, epsrel=epsrel) s *= 6.2831853071795862 ds *= 6.2831853071795862 return s, ds if __name__ == '__main__': ftel = 22.5 # Telescope focal length [m] pixsize = 15e-6 # Pixel size [m] px2arcs = pixsize/ftel*RAD2ARCS # Pixel scale [arcsec/px] print "Focal length=%.2f m, pixel size=%.0f microns, scale=%.3f\"/px" % \ (ftel, pixsize*1e6, px2arcs) amax = 20. # Image half width [arcsec] xmax = int(amax / px2arcs) # Image half width [px] amax = xmax * px2arcs # Image half width [arcsec] print "Image half width: %d px = %.2f\"" % (xmax, amax) x = N.arange(-xmax, xmax+0.1, dtype='d') # (n,) [px] r = N.hypot(*N.meshgrid(x, x)) # Radius 2D-array (n,n) [px] a = r * px2arcs # Radius (n,n) [arcsec] px2f = ftel/pixsize/len(x) # [px/rad] f = r * px2f # Spatial freq 2D-array (n,n) [1/rad] # Reference medium seeing lref = 5e-7 # Reference wavelength [m] r0ref = 0.10 # Fried parameter at ref. lbda [m] expo = 5/3 # expo = 1.5 # Kolmogorov requires 5/3 # Medium seeing at 5000A lbda = 5e-7 # Current wavelength [m] r0m = friedParamater(lbda, r0ref, lref, expo) s = seeing_fwhm(lbda, r0ref, lref, expo, verbose=True) # Seeing FWHM [arcsec] psf = psf2D_FFT(otf_Kolmogorov(lbda*f, r0=r0m, expo=expo)) ar = a.ravel() # (n**2,) j = N.argsort(ar) i = j[:N.searchsorted(ar[j], 5)] # ar[i]: r up to 5 arcsec r1 = N.logspace(-1, N.log10(5/px2arcs), 100) # Radius (n,) [px] a1 = r1 * px2arcs # Radius (n,) [arcsec] psf1 = psf_Kolmogorov_Hankel(a1, lbda, r0=r0m, expo=expo) psf2 = psf_Kolmogorov_FFT(a1, lbda, scale=px2arcs, r0=r0m, expo=expo) print "Integral(Laguerre)", \ integ_radial_laguerre( lambda r: psf_Kolmogorov_Hankel(r, lbda, r0=r0m, expo=expo)) print "Integral(quad)", \ integ_radial_quad( lambda r: psf_Kolmogorov_Hankel(r, lbda, r0=r0m, expo=expo)) if False: fig = P.figure() ax = fig.add_subplot(1, 1, 1, title="Kolmogorov PSF (expo=%.2f)" % expo, xlabel="r [arcsec]", xscale='log', ylabel="Normalized flux", yscale='log', ) ax.plot(ar[i], psf.ravel()[i], color=RED, ls='None', marker='.', ms=5, label='2D-FFT') ax.plot(a1, psf2, color=GREEN, ls='None', marker='.', ms=5, label='Interpolated 2D-FFT') ax.plot(a1, psf1, label='Hankel') #ax.plot(a1,a1*psf1*N.exp(a1), label='Hankel') ax.axhline(0.5, c='0.8', label='_') ax.axvline(s/2, c='0.8', label='_') ax.set_ylim(psf1.min()/1.1, psf1.max()*1.1) ax.legend(loc='best', fontsize='small') # Plots ============================== if False: fig = P.figure() axI = fig.add_subplot(1, 1, 1, title=u"PSF (r₀=%.0f cm@%.0f Å)" % (r0m*100, lbda*1e10), xlabel='x [arcsec]', ylabel='y [arcsec]', aspect='equal') psf = psf_Kolmogorov_Hankel(a, lbda, r0=r0m) print psf.sum()/len(a)**2 print psf.max(), lbda, S.gamma(6/5) axI.imshow(psf, norm=P.matplotlib.colors.LogNorm(vmin=(psf[psf > 0]).min()), interpolation='nearest', extent=[-amax, amax, -amax, amax]) # Radial profile ------------------------------ if True: fig = P.figure() axR = fig.add_subplot(1, 1, 1, title=u"Kolmogorov seeing radial profile " u"(r₀=%.0f cm@%.0f Å)" % (r0ref*1e2, lref*1e10), xlabel='r [arcsec]', # xscale='log', ylabel='Normalized flux', yscale='log') axR.plot(ar[i], psf.ravel()[i], color=GREEN, ls='-', marker='None', label=u'%.0f Å: seeing=%.2f"' % (lbda*1e10, s)) lbda = 3e-7 # Medium seeing @3000A r0 = friedParamater(lbda, r0ref, lref) s = seeing_fwhm(lbda, r0ref, lref) axR.plot(ar[i], psf2D_FFT(otf_Kolmogorov(lbda*f, r0)).ravel()[i], color=BLUE, ls='-', marker='None', label=u'%.0f Å: seeing=%.2f"' % (lbda*1e10, s)) lbda = 10e-7 # Medium seeing @10000A r0 = friedParamater(lbda, r0ref, lref) s = seeing_fwhm(lbda, r0ref, lref) axR.plot(ar[i], psf2D_FFT(otf_Kolmogorov(lbda*f, r0)).ravel()[i], color=RED, ls='-', marker='None', label=u'%.0f Å: seeing=%.2f"' % (lbda*1e10, s)) lbda = 5e-7 r0ref = 0.15 # Good seeing @5000A r0 = friedParamater(lbda, r0ref, lref) s = seeing_fwhm(lbda, r0ref, lref) axR.plot(ar[i], psf2D_FFT(otf_Kolmogorov(lbda*f, r0)).ravel()[i], color=GREEN, ls=':', marker='None', label=u'r₀=%.0f cm: seeing=%.2f"' % (r0ref*100, s)) r0ref = 0.07 # Bad seeing @5000A r0 = friedParamater(lbda, r0ref, lref) s = seeing_fwhm(lbda, r0ref, lref) axR.plot(ar[i], psf2D_FFT(otf_Kolmogorov(lbda*f, r0)).ravel()[i], color=GREEN, ls='--', marker='None', label=u'r₀=%.0f cm: seeing=%.2f"' % (r0ref*100, s)) axR.set_ylim(1e-4, 1) axR.legend(loc='best', fontsize='small') # Comparison w/ extract_star PSF ------------------------------ if True: lbda = 5e-7 r0ref = 0.10 s = seeing_fwhm(lbda, r0ref, lref) fig = P.figure() ax = fig.add_subplot(1, 1, 1, title=u"Kolmogorov (r₀=%.0f cm@%.0f Å) " "vs. ad-hoc seeing PSF" % (r0ref*1e2, lref*1e10), xlabel='r [arcsec]', # xscale='log', ylabel='Normalized flux', yscale='log') # Long exposures gl, ml = psf_ES_long(ar[i], 2.03, decompose=True) ax.plot(ar[i], gl+ml, color=BLUE, ls='-', marker='None', label="Ad-hoc PSF (long)") ax.set_autoscale_on(False) ax.plot(ar[i], gl, color=BLUE, ls='--', marker='None', label="_") ax.plot(ar[i], ml, color=BLUE, ls='--', marker='None', label="_") # Short exposures gs, ms = psf_ES_short(ar[i], 2, decompose=True) ax.plot(ar[i], gs+ms, color=RED, ls='-', marker='None', label="Ad-hoc PSF (short)") ax.plot(ar[i], gs, color=RED, ls='--', marker='None', label="_") # [Individual components] ax.plot(ar[i], ms, color=RED, ls='--', marker='None', label="_") # Kolmogorov ax.plot(ar[i], psf.ravel()[i], color=GREEN, ls='-', lw=2, marker='None', label=u'%.0f Å: seeing=%.2f"' % (lbda*1e10, s)) for n in range(8, 12): y = psf_Kolmogorov_Hankel(a1, lbda, r0=r0m, expo=n/6) if n != 10: ax.plot(a1[72:], y[72:], color=GREEN, ls='-', marker='None', label='_') ax.annotate('n=%s' % (fractions.Fraction(n, 6)), (a1[90], y[90]), fontsize='small', rotation=-25) ax.set_ylim(1e-4, 1) ax.legend(loc='best', fontsize='small') # Ad-hoc comparison w/ observations ------------------------------ if True: # Long exposure @4330A: GD71_07_253_090_B.txt Rl, Fl = N.loadtxt("GD71_07_253_090_B.txt", unpack=True) Fl /= Fl.max()/0.8 # /0.85 ll = 4.33e-7 rl = 0.038 psfl = psf2D_FFT(otf_Kolmogorov(ll*f, rl)) + 0.0035 psfl2 = psf2D_FFT(otf_Kolmogorov(ll*f, rl, expo=1.52)) + 0.0035 # Short exposure @7600: HR7950_07_253_041_R.txt Rs, Fs = N.loadtxt("HR7950_07_253_041_R.txt", unpack=True) Fs /= Fs.max()/0.9 ls = 7.6e-7 rs = 0.075 psfs = psf2D_FFT(otf_Kolmogorov(ls*f, rs)) psfs2 = psf2D_FFT(otf_Kolmogorov(ls*f, rs, expo=1.5)) fig = P.figure(figsize=(10, 5)) axl = fig.add_subplot(1, 2, 1, title=u"Long exposure @%.0f Å" % (ll*1e10), xlabel="r [arcsec]", ylabel="Normalized flux", yscale='log') axs = fig.add_subplot(1, 2, 2, title=u"Short exposure @%.0f Å" % (ls*1e10), xlabel="r [arcsec]", yscale='log') i = j[:N.searchsorted(ar[j], 10)] # Keep r up to 10 arcsec axl.plot(Rl, Fl, color=BLUE, ls='None', marker='.', ms=5, label='Observed') axl.plot(ar[i], psfl.ravel()[i], color=GREEN, ls='-', lw=2, marker='None', label='Kolmogorov (n=5/3)', zorder=0) axl.plot(ar[i], psfl2.ravel()[i], color=GREEN, ls='-', marker='None', label='Kolmogorov (n=3/2)', zorder=0) axs.plot(Rs, Fs, color=RED, ls='None', marker='.', ms=5, label='Observed') axs.plot(ar[i], psfs.ravel()[i], color=GREEN, ls='-', lw=2, marker='None', label='Kolmogorov (n=5/3)', zorder=0) axs.plot(ar[i], psfs2.ravel()[i], color=GREEN, ls='-', marker='None', label='Kolmogorov (n=3/2)', zorder=0) axl.set_ylim(1e-4, 1.1) axs.set_ylim(1e-4, 1.1) axl.legend(loc='best', fontsize='small') axs.legend(loc='best', fontsize='small') if True: for i in P.get_fignums(): fig = P.figure(i) fig.savefig("seeing_%d.png" % i) P.show()