from __future__ import division import sys import os import getopt import math import random import re, glob from math import sin, cos, floor, sqrt, pi, tan, atan # asin, atan2 from scipy.stats import distributions from scipy.stats import chisquare from mx.DateTime import * import bisect try: import DateTime except ImportError: from mx import DateTime #******************** Moon tools ******************* # Precision used when describing the moon's phase in textual format, # in phase_string(). PRECISION = 0.05 NEW = 0 / 4.0 FIRST = 1 / 4.0 FULL = 2 / 4.0 LAST = 3 / 4.0 NEXTNEW = 4 / 4.0 class MoonPhase: """I describe the phase of the moon. I have the following properties: date - a DateTime instance phase - my phase, in the range 0.0 .. 1.0 phase_text - a string describing my phase illuminated - the percentage of the face of the moon illuminated angular_diameter - as seen from Earth, in degrees. sun_angular_diameter - as seen from Earth, in degrees. new_date - the date of the most recent new moon q1_date - the date the moon reaches 1st quarter in this cycle full_date - the date of the full moon in this cycle q3_date - the date the moon reaches 3rd quarter in this cycle nextnew_date - the date of the next new moon """ def __init__(self, date=DateTime.now()): """MoonPhase constructor. Give me a date, as either a Julian Day Number or a DateTime object.""" if not isinstance(date, DateTime.DateTimeType): self.date = DateTime.DateTimeFromJDN(date) else: self.date = date self.__dict__.update(phase(self.date)) self.phase_text = phase_string(self.phase) def __getattr__(self, a): if a in ['new_date', 'q1_date', 'full_date', 'q3_date', 'nextnew_date']: (self.new_date, self.q1_date, self.full_date, self.q3_date, self.nextnew_date) = phase_hunt(self.date) return getattr(self,a) raise AttributeError(a) def __repr__(self): if type(self.date) is int: jdn = self.date else: jdn = self.date.jdn return "<%s(%d)>" % (self.__class__, jdn) def __str__(self): if type(self.date) is int: d = DateTime.DateTimeFromJDN(self.date) else: d = self.date s = "%s for %s, %s (%%%.2f illuminated)" %\ (self.__class__, d.strftime(), self.phase_text, self.illuminated * 100) return s class AstronomicalConstants: # JDN stands for Julian Day Number # Angles here are in degrees # 1980 January 0.0 in JDN # XXX: DateTime(1980).jdn yields 2444239.5 -- which one is right? epoch = 2444238.5 # Ecliptic longitude of the Sun at epoch 1980.0 ecliptic_longitude_epoch = 278.833540 # Ecliptic longitude of the Sun at perigee ecliptic_longitude_perigee = 282.596403 # Eccentricity of Earth's orbit eccentricity = 0.016718 # Semi-major axis of Earth's orbit, in kilometers sun_smaxis = 1.49585e8 # Sun's angular size, in degrees, at semi-major axis distance sun_angular_size_smaxis = 0.533128 ## Elements of the Moon's orbit, epoch 1980.0 # Moon's mean longitude at the epoch moon_mean_longitude_epoch = 64.975464 # Mean longitude of the perigee at the epoch moon_mean_perigee_epoch = 349.383063 # Mean longitude of the node at the epoch node_mean_longitude_epoch = 151.950429 # Inclination of the Moon's orbit moon_inclination = 5.145396 # Eccentricity of the Moon's orbit moon_eccentricity = 0.054900 # Moon's angular size at distance a from Earth moon_angular_size = 0.5181 # Semi-mojor axis of the Moon's orbit, in kilometers moon_smaxis = 384401.0 # Parallax at a distance a from Earth moon_parallax = 0.9507 # Synodic month (new Moon to new Moon), in days synodic_month = 29.53058868 # Base date for E. W. Brown's numbered series of lunations (1923 January 16) lunations_base = 2423436.0 ## Properties of the Earth earth_radius = 6378.16 c = AstronomicalConstants() # Little handy mathematical functions. fixangle = lambda a: a - 360.0 * floor(a/360.0) torad = lambda d: d * pi / 180.0 todeg = lambda r: r * 180.0 / pi dsin = lambda d: sin(torad(d)) dcos = lambda d: cos(torad(d)) def phase_string(p): phase_strings = ( (NEW + PRECISION, "new"), (FIRST - PRECISION, "waxing crescent"), (FIRST + PRECISION, "first quarter"), (FULL - PRECISION, "waxing gibbous"), (FULL + PRECISION, "full"), (LAST - PRECISION, "waning gibbous"), (LAST + PRECISION, "last quarter"), (NEXTNEW - PRECISION, "waning crescent"), (NEXTNEW + PRECISION, "new")) i = bisect.bisect([a[0] for a in phase_strings], p) return phase_strings[i][1] def phase(phase_date=DateTime.now()): """Calculate phase of moon as a fraction: The argument is the time for which the phase is requested, expressed in either a DateTime or by Julian Day Number. Returns a dictionary containing the terminator phase angle as a percentage of a full circle (i.e., 0 to 1), the illuminated fraction of the Moon's disc, the Moon's age in days and fraction, the distance of the Moon from the centre of the Earth, and the angular diameter subtended by the Moon as seen by an observer at the centre of the Earth.""" # Calculation of the Sun's position # date within the epoch if hasattr(phase_date, "jdn"): day = phase_date.jdn - c.epoch else: day = phase_date - c.epoch # Mean anomaly of the Sun N = fixangle((360/365.2422) * day) # Convert from perigee coordinates to epoch 1980 M = fixangle(N + c.ecliptic_longitude_epoch - c.ecliptic_longitude_perigee) # Solve Kepler's equation Ec = kepler(M, c.eccentricity) Ec = sqrt((1 + c.eccentricity) / (1 - c.eccentricity)) * tan(Ec/2.0) # True anomaly Ec = 2 * todeg(atan(Ec)) # Suns's geometric ecliptic longuitude lambda_sun = fixangle(Ec + c.ecliptic_longitude_perigee) # Orbital distance factor F = ((1 + c.eccentricity * cos(torad(Ec))) / (1 - c.eccentricity**2)) # Distance to Sun in km sun_dist = c.sun_smaxis / F sun_angular_diameter = F * c.sun_angular_size_smaxis ######## # # Calculation of the Moon's position # Moon's mean longitude moon_longitude = fixangle(13.1763966 * day + c.moon_mean_longitude_epoch) # Moon's mean anomaly MM = fixangle(moon_longitude - 0.1114041 * day - c.moon_mean_perigee_epoch) # Moon's ascending node mean longitude # MN = fixangle(c.node_mean_longitude_epoch - 0.0529539 * day) evection = 1.2739 * sin(torad(2*(moon_longitude - lambda_sun) - MM)) # Annual equation annual_eq = 0.1858 * sin(torad(M)) # Correction term A3 = 0.37 * sin(torad(M)) MmP = MM + evection - annual_eq - A3 # Correction for the equation of the centre mEc = 6.2886 * sin(torad(MmP)) # Another correction term A4 = 0.214 * sin(torad(2 * MmP)) # Corrected longitude lP = moon_longitude + evection + mEc - annual_eq + A4 # Variation variation = 0.6583 * sin(torad(2*(lP - lambda_sun))) # True longitude lPP = lP + variation # # Calculation of the Moon's inclination # unused for phase calculation. # Corrected longitude of the node # NP = MN - 0.16 * sin(torad(M)) # Y inclination coordinate # y = sin(torad(lPP - NP)) * cos(torad(c.moon_inclination)) # X inclination coordinate # x = cos(torad(lPP - NP)) # Ecliptic longitude (unused?) # lambda_moon = todeg(atan2(y,x)) + NP # Ecliptic latitude (unused?) # BetaM = todeg(asin(sin(torad(lPP - NP)) * sin(torad(c.moon_inclination)))) ####### # # Calculation of the phase of the Moon # Age of the Moon, in degrees moon_age = lPP - lambda_sun # Phase of the Moon moon_phase = (1 - cos(torad(moon_age))) / 2.0 # Calculate distance of Moon from the centre of the Earth moon_dist = (c.moon_smaxis * (1 - c.moon_eccentricity**2))\ / (1 + c.moon_eccentricity * cos(torad(MmP + mEc))) # Calculate Moon's angular diameter moon_diam_frac = moon_dist / c.moon_smaxis moon_angular_diameter = c.moon_angular_size / moon_diam_frac # Calculate Moon's parallax (unused?) # moon_parallax = c.moon_parallax / moon_diam_frac res = { 'phase': fixangle(moon_age) / 360.0, 'illuminated': moon_phase, 'age': c.synodic_month * fixangle(moon_age) / 360.0 , 'distance': moon_dist, 'angular_diameter': moon_angular_diameter, 'sun_distance': sun_dist, 'sun_angular_diameter': sun_angular_diameter } return res # phase() def phase_hunt(sdate=DateTime.now()): """Find time of phases of the moon which surround the current date. Five phases are found, starting and ending with the new moons which bound the current lunation. """ if not hasattr(sdate,'jdn'): sdate = DateTime.DateTimeFromJDN(sdate) adate = sdate + DateTime.RelativeDateTime(days=-45) k1 = floor((adate.year + ((adate.month - 1) * (1.0/12.0)) - 1900) * 12.3685) nt1 = meanphase(adate, k1) adate = nt1 sdate = sdate.jdn while 1: adate = adate + c.synodic_month k2 = k1 + 1 nt2 = meanphase(adate,k2) if nt1 <= sdate < nt2: break nt1 = nt2 k1 = k2 phases = list(map(truephase, [k1, k1, k1, k1, k2], [0/4.0, 1/4.0, 2/4.0, 3/4.0, 0/4.0])) return phases # phase_hunt() def meanphase(sdate, k): """Calculates time of the mean new Moon for a given base date. This argument K to this function is the precomputed synodic month index, given by: K = (year - 1900) * 12.3685 where year is expressed as a year and fractional year. """ # Time in Julian centuries from 1900 January 0.5 if not hasattr(sdate,'jdn'): delta_t = sdate - DateTime.DateTime(1900,1,1,12).jdn t = delta_t / 36525 else: delta_t = sdate - DateTime.DateTime(1900,1,1,12) t = delta_t.days / 36525 # square for frequent use t2 = t * t # and cube t3 = t2 * t nt1 = ( 2415020.75933 + c.synodic_month * k + 0.0001178 * t2 - 0.000000155 * t3 + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2) ) return nt1 # meanphase() def truephase(k, tphase): """Given a K value used to determine the mean phase of the new moon, and a phase selector (0.0, 0.25, 0.5, 0.75), obtain the true, corrected phase time.""" apcor = False # add phase to new moon time k = k + tphase # Time in Julian centuries from 1900 January 0.5 t = k / 1236.85 t2 = t * t t3 = t2 * t # Mean time of phase pt = ( 2415020.75933 + c.synodic_month * k + 0.0001178 * t2 - 0.000000155 * t3 + 0.00033 * dsin(166.56 + 132.87 * t - 0.009173 * t2) ) # Sun's mean anomaly m = 359.2242 + 29.10535608 * k - 0.0000333 * t2 - 0.00000347 * t3 # Moon's mean anomaly mprime = 306.0253 + 385.81691806 * k + 0.0107306 * t2 + 0.00001236 * t3 # Moon's argument of latitude f = 21.2964 + 390.67050646 * k - 0.0016528 * t2 - 0.00000239 * t3 if (tphase < 0.01) or (abs(tphase - 0.5) < 0.01): # Corrections for New and Full Moon pt = pt + ( (0.1734 - 0.000393 * t) * dsin(m) + 0.0021 * dsin(2 * m) - 0.4068 * dsin(mprime) + 0.0161 * dsin(2 * mprime) - 0.0004 * dsin(3 * mprime) + 0.0104 * dsin(2 * f) - 0.0051 * dsin(m + mprime) - 0.0074 * dsin(m - mprime) + 0.0004 * dsin(2 * f + m) - 0.0004 * dsin(2 * f - m) - 0.0006 * dsin(2 * f + mprime) + 0.0010 * dsin(2 * f - mprime) + 0.0005 * dsin(m + 2 * mprime) ) apcor = True elif (abs(tphase - 0.25) < 0.01) or (abs(tphase - 0.75) < 0.01): pt = pt + ( (0.1721 - 0.0004 * t) * dsin(m) + 0.0021 * dsin(2 * m) - 0.6280 * dsin(mprime) + 0.0089 * dsin(2 * mprime) - 0.0004 * dsin(3 * mprime) + 0.0079 * dsin(2 * f) - 0.0119 * dsin(m + mprime) - 0.0047 * dsin(m - mprime) + 0.0003 * dsin(2 * f + m) - 0.0004 * dsin(2 * f - m) - 0.0006 * dsin(2 * f + mprime) + 0.0021 * dsin(2 * f - mprime) + 0.0003 * dsin(m + 2 * mprime) + 0.0004 * dsin(m - 2 * mprime) - 0.0003 * dsin(2 * m + mprime) ) if (tphase < 0.5): # First quarter correction pt = pt + 0.0028 - 0.0004 * dcos(m) + 0.0003 * dcos(mprime) else: # Last quarter correction pt = pt + -0.0028 + 0.0004 * dcos(m) - 0.0003 * dcos(mprime) apcor = True if not apcor: raise ValueError( "TRUEPHASE called with invalid phase selector", tphase) return DateTime.DateTimeFromJDN(pt) # truephase() def kepler(m, ecc): """Solve the equation of Kepler.""" epsilon = 1e-6 m = torad(m) e = m while 1: delta = e - ecc * sin(e) - m e = e - delta / (1.0 - ecc * cos(e)) if abs(delta) <= epsilon: break return e #******************** chi square ******************* def GammaInc_Q( a, x): a1 = a-1 a2 = a-2 def f0( t ): return t**a1*math.exp(-t) def df0(t): return (a1-t)*t**a2*math.exp(-t) y = a1 while f0(y)*(x-y) >2.0e-8 and y < x: y += .3 if y > x: y = x h = 3.0e-4 n = int(y/h) h = y/n hh = 0.5*h gamax = h * sum( f0(t)+hh*df0(t) for t in ( h*j for j in xrange(n-1, -1, -1))) return gamax/gamma_spounge(a) co = None def gamma_spounge( z): global co a = 12 if co is None: k1_factrl = 1.0 co = [] co.append(math.sqrt(2.0*math.pi)) for k in range(1,a): co.append( math.exp(a-k) * (a-k)**(k-0.5) / k1_factrl ) k1_factrl *= -k accm = co[0] for k in range(1,a): accm += co[k] / (z+k) accm *= math.exp( -(z+a)) * (z+a)**(z+0.5) return accm/z; def chi2UniformDistance( dataSet ): expected = sum(dataSet)*1.0/len(dataSet) cntrd = (d-expected for d in dataSet) return sum(x*x for x in cntrd)/expected def chi2Probability(dof, distance): return 1.0 - GammaInc_Q( 0.5*dof, 0.5*distance) def chi2IsUniform(dataSet, significance): dof = len(dataSet)-1 dist = chi2UniformDistance(dataSet) return chi2Probability( dof, dist ) > significance #################################################################### def usage(): print """ This is the help docstring of DateParserCIP. This is a (very simple and sometimes even candid!) Python script that i) extracts dates in ISO format (yyyy-mm-dd) from a (series of) input textual file(s) whose name must be provided as argument from the command line with option -i, or --input, and that ii) computes how these dates are distributed within the solar year, the week and the sinodic month (i.e., moon phases). Input options: -h :: to display this help -i :: to specify the input file -s :: to specify the first date of the time interval to consider (e.g., -s 1974-04-18) -u :: to specify the last date of the time interval to consider (e.g., -u 2012-12-21) Output: The script produces four txt files (and puts them in the root directory): - output.txt :: it is the main output file with high level statistics and some aggregate results. - PleniNoviList.txt :: it is a sort of check file, as it reports the novilunes and plenilunes found in the considered interval, so that they can be checked with e.g., http://eclipse.gsfc.nasa.gov/phase/phases1901.html - serietemporale.txt :: it is the time series of the number of event occurrences in each specific date. It can be huge. - matricesinodica.txt :: it is the horizontal and vertical representation of the number of occurrences for each single date along the lunar month. It can be even huger. Release Notes: - Mind that the time interval (if specified) must be of at least 100 days (otherwise, the progress indicator gives a... typical trouble) - This script uses a very accurate algorithm for moon age calculation, ported in python by Kevin Turner (http://bazaar.launchpad.net/~keturn/py-moon-phase/trunk/annotate/head:/moon.py) on the basis of a moontool by John Walker (http://www.fourmilab.ch/). For this reason, this script can be used to assess if any moon effect is detectable in the above mentioned distributions. Indeed, I wrote this script right to detect the moon effect on birth dates. - Mind that I did not bother with checking if inputs are well formatted and if they even exist; so be careful (re exception raisings) and accurate (re feedings)! - Also mind that I did not bother with performance optmization (should I have?); so be understanding (re me) and patient (re the processing time)! - This script is released 'as-is' under the Creative Commons "license" CC BY-NC-SA http://creativecommons.org/licenses/by-nc-sa/3.0/. Enjoy! FC 2011-03-07, Milan (Italy) """ def main(): input = '' datesincestringa = '' dateuntilstringa = '' try: options, remainder = getopt.getopt(sys.argv[1:], 'hi:s:u:', ["input=", "help", "since", "until"]) except getopt.GetoptError, err: # print help information and exit: print str(err) # will print something like "option -a not recognized" usage() sys.exit(2) for o, a in options: #print 'o : %s; a: %s' % (o, a) if o in ("-i", "--input"): input = a elif o in ("-h", "--help"): usage() sys.exit() elif o in ("-s", "--since"): datesincestringa = a elif o in ("-u", "--until"): dateuntilstringa = a else: assert False, "unhandled option" # ****************** COSTANTI INIZIALI *************************** distribuzioneAssoluta = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] distribuzioneProporzioni = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] distribuzioneGestazione = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] distribuzionemesi = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] fasi = {'new':0, 'waxing crescent':0, 'first quarter':0, 'waxing gibbous':0, 'full':0, 'waning gibbous':0, 'last quarter':0, 'waning crescent':0, 'new':0} week = [0, 0, 0, 0, 0, 0, 0] dateveramenteconsiderate = [0] sinodico = [[0] for i in range(30)] sinodicofasi = {'new':[0], 'waxing crescent':[0], 'first quarter':[0], 'waxing gibbous':[0], 'full':[0], 'waning gibbous':[0], 'last quarter':[0], 'waning crescent':[0], 'new':[0]} perigees = [[]] apogees = [[]] path = './%s' % (input) # ****************** ELABORAZIONE *************************** inputfile = '' if path == './': #gestione scanning di tutti i file txt nella directory indicefile = 0 contenutofiles = '' for infile in glob.glob( os.path.join(path, '*.txt') ): indicefile = indicefile + 1 print str(indicefile) + ") Processing file: " + infile inputfile = inputfile + str(infile) + '\n' try: data = open(infile,'r') contenutofile = data.read() except: contenutofile = '' contenutofiles = contenutofiles + contenutofile data.close() datedinascita = re.findall(r'\d{4}-\d{2}-\d{2}', contenutofiles) else: inputfile = input try: data = open(path,'r') contenutofile = data.read() datedinascita = re.findall(r'\d{4}-\d{2}-\d{2}', contenutofile) settimane = re.findall(r',\w{2}', contenutofile) sessi = re.findall(r',\w{1},', contenutofile) data.close() except: print "Error in opening the input file." sys.exit() noduplicates = list(set(datedinascita)) numerodate = len(noduplicates) serietemporale=[0] datesbagliate = 0 numeromaschi = 0 anno, mese, giorno = datedinascita[0].split('-') try: dataminima = DateTime.Date(int(anno), int(mese), int(giorno)) except: dataminima = DateTime.Date(1974, 1, 1) try: datamassima = DateTime.Date(int(anno), int(mese), int(giorno)) except: datamassima = DateTime.Date(1974, 1, 1) if datesincestringa == '': datesince = DateTime.Date(1880, 1, 1) else: annosince, mesesince, giornosince = datesincestringa.split('-') datesince = DateTime.Date(int(annosince), int(mesesince), int(giornosince)) if dateuntilstringa == '': dateuntil = DateTime.now() else: annountil, meseuntil, giornountil = dateuntilstringa.split('-') dateuntil = DateTime.Date(int(annountil), int(meseuntil), int(giornountil)) inizio = DateTime.now() numeroeventi = len(datedinascita) for i in range(0, numeroeventi): year, month, day = datedinascita[i].split('-') try: datanascita = DateTime.Date(int(year), int(month), int(day), 10, 30) #date is considered at 12:30 am UMT + 1 to make a daily average and consider daylight saving times (sort of). if ((datanascita > DateTime.Date(1861, 3, 17)) and (datanascita < DateTime.now())): if datanascita <= dataminima: dataminima = datanascita if datanascita >= datamassima: datamassima = datanascita if (datanascita >= datesince) and (datanascita <= dateuntil): #numeronatioggi = datedinascita.count(datedinascita[i]) dateveramenteconsiderate.append(datedinascita[i]) m = MoonPhase(datanascita) ageofthemoon = m.age giornolunare = int(ageofthemoon) giornosettimana = datanascita.weekday() #print '%d) per la data %s associata eta` lunare %d' % (i, datanascita, giornolunare) distribuzioneAssoluta[giornolunare] +=1 fasi[m.phase_text] +=1 week[giornosettimana] +=1 distribuzionemesi[int(month)] +=1 else: datesbagliate += 1 if (i % (numeroeventi // 100)) == 0: print 'Sinodic Distribution :: Done %.2f percent' % (((i+1)/numeroeventi)*100) # except ValueError: datesbagliate += 1 pass print "Do not Panic, date %s is not valid" % (datedinascita[i]) if dateuntilstringa == '': dateuntil = datamassima if datesincestringa == '': datesince = dataminima dateveramenteconsiderate.remove(0) numeronascite = len(dateveramenteconsiderate) noduplicatesininterval = list(set(dateveramenteconsiderate)) #if (dateuntilstringa <> '') and (datesincestringa <> ''): numerodateconsiderate = (len(noduplicatesininterval)) #else: # numerodateconsiderate = numerodate if len(settimane) > 2: settimanepresenti = True for i in range(0, len(settimane)): settimana = int(settimane[i].strip(',')) distribuzioneGestazione[settimana] +=1 else: settimanepresenti = False for i in range(0, len(sessi)): sesso = sessi[i].strip(',') if sesso == 'M': numeromaschi +=1 secondiintervallo = int(datamassima - dataminima) giorniintervallo = (secondiintervallo / 86400) + 1 secondiintervalloconsiderato = int(dateuntil - datesince) giorniintervalloconsiderato = (secondiintervalloconsiderato / 86400) + 1 prog = 0 elencopleni = '' elenconovi = '' for elem in noduplicatesininterval: numdate = dateveramenteconsiderate.count(elem) #if numdate > 100: # print "OUTLIER ----> %s with %d events" % (elem, numdate) anno, mese, giorno = elem.split('-') try: datanascita = DateTime.Date(int(anno), int(mese), int(giorno), 10, 30) except ValueError: pass luna = MoonPhase(datanascita) giornolunare = int(luna.age) sinodico[giornolunare].append(numdate) fase = luna.phase_text sinodicofasi[fase].append(numdate) distanza = luna.distance if distanza < 365000: perigees.append([elem,numdate]) if distanza > 405000: apogees.append([elem,numdate]) #sinodico[giornolunare][indice[giornolunare]] = numdate #indice[giornolunare] +=1 prog += 1 if fase == 'full': elencopleni = elencopleni + str(datanascita) + ', ' + str(giornolunare) + '\n' if fase == 'new': elenconovi = elenconovi + str(datanascita) + ', ' + str(giornolunare) + '\n' if (prog % (numerodateconsiderate // 100)) == 0: print "Sinodic Matrix :: Done %.2f percent" % ((prog / numerodateconsiderate)*100) fine = DateTime.now() durata = (fine - inizio) # # ****************** SCRITTURA OUTPUT *************************** timestamp = fine.strftime('%Y%m%d%H%M') print 'Now writing the output file' path = './output' + timestamp + '.txt' dataoutput = open(path,'w') dataoutput.write( "File(s) processed: " + inputfile + '\n') dataoutput.write( 'Processing time: %.2f seconds. \n' % (durata) ) #print 'numerodate: %d, numerodateconsiderate: %d, giorniintervalloconsiderato: %d, numeronascite: %d' % (numerodate, numerodateconsiderate, giorniintervalloconsiderato, numeronascite) dataoutput.write( "Total number of dates found in input file: %d within the time period since %s until %s. \n" % (numerodate, dataminima.strftime('%Y-%m-%d'), datamassima.strftime('%Y-%m-%d')) ) dataoutput.write( "Total number of dates processed: %d (on %d potential dates within the time period " % ((numerodateconsiderate + 1), giorniintervalloconsiderato) ) if datesincestringa <> '': dataoutput.write( "since %s " % (datesince.strftime('%Y-%m-%d')) ) if dateuntilstringa <> '': dataoutput.write( "until %s.)" % (dateuntil.strftime('%Y-%m-%d')) ) dataoutput.write( '\n Total number of events (e.g., births): %d, on average %.2f on each day. \n\n' % (numeronascite, (numeronascite / numerodateconsiderate)) ) iq = ((datesbagliate / numeronascite) * 100) if iq <> 0: dataoutput.write( 'IQ level (i.e., percentage of dates wrongly formatted): %.3f \n\n' % ( iq )) if numeromaschi <> 0: dataoutput.write( 'Male population proportion: %.2f \n\n' % ((numeromaschi / numeronascite) * 100) ) dataoutput.write( 'Distribution of events within the sinodic month: \n' ) for el in range(0, len(distribuzioneProporzioni)): prop = (distribuzioneAssoluta[el] / numeronascite) distribuzioneProporzioni[el] = prop * 100 prod = (1-prop) / numeronascite sep = sqrt(prop * prod) #vone = 2 * distribuzioneAssoluta[el] #vtwo = 2 * (numeronascite - distribuzioneAssoluta[el] + 1) #vthree = 2 * (distribuzioneAssoluta[el] + 1) #vfour = 2 * (numeronascite - distribuzioneAssoluta[el]) #u = distributions.f.ppf(0.95, vone, vtwo) #w = distributions.f.ppf(0.95, vthree, vfour) #lower = (vone * u) / (vtwo + (vone * u)) #upper = (vthree * w) / (vfour + (vthree * w)) dataoutput.write( 'Day: %d, event no.: %d, percentage: %.2f, lowerbound: %.2f, upperbound: %.2f \n' % (el, distribuzioneAssoluta[el], distribuzioneProporzioni[el], (prop - (1.96 * sep))*100, (prop + (1.96 * sep))*100) ) #dataoutput.write( 'Day: %d, event no.: %d, percentage: %.2f, lb: %.2f, ub: %.2f, lowerbound: %.2f, upperbound: %.2f \n' % (el, distribuzioneAssoluta[el], distribuzioneProporzioni[el], (prop - (1.96 * sep))*100, (prop + (1.96 * sep))*100, ((prop - lower) * 100), ((prop + upper) * 100)) ) #for el in range(0, len(distribuzioneProporzioni)): # distribuzioneProporzioni[el] = ((distribuzioneAssoluta[el] / numeronascite) * 100) # dataoutput.write( 'Day: %d, event no.: %d, percentage: %.2f \n' % (el, distribuzioneAssoluta[el], distribuzioneProporzioni[el]) ) dof = len(distribuzioneAssoluta)-1 distance = chi2UniformDistance(distribuzioneAssoluta) prob = chi2Probability( dof, distance) dataoutput.write('') dataoutput.write( "Chi Square: dof: %d, distance: %.4f, P-value: %.4f \n\n" % (dof, distance, prob) ) dataoutput.write( 'Distribution of events within the week (MON-SUN): \n' ) for i in range(0, len(week)): dataoutput.write( 'Day %d: %d (%.2f)\n' % (i, week[i], ((week[i]/numeronascite)*100)) ) dataoutput.write( '' ) dof = len(week)-1 distance =chi2UniformDistance(week) prob = chi2Probability( dof, distance) dataoutput.write( "Chi Square: dof: %d, distance: %.4f, P-value: %.4f \n\n" % (dof, distance, prob) ) dataoutput.write( "Distribution of events within the solar year (JAN-DEC):\n" ) for i in range(1, len(distribuzionemesi)): dataoutput.write( 'Month %d: %d (%.2f)\n' % (i, distribuzionemesi[i], ((distribuzionemesi[i]/numeronascite)*100)) ) dataoutput.write( '' ) dof = len(distribuzionemesi)-1 distance =chi2UniformDistance(distribuzionemesi) prob = chi2Probability( dof, distance) dataoutput.write( "Chi Square: dof: %d, distance: %.4f, P-value: %.4f \n\n" % (dof, distance, prob) ) if settimanepresenti: dataoutput.write( 'Distribution of events by Gestation Week (since the 20th on): \n' ) for i in range(20, len(distribuzioneGestazione)): dataoutput.write( 'Week %d: %d, (%.2f)\n' % (i, distribuzioneGestazione[i], ((distribuzioneGestazione[i]/numeronascite)*100)) ) dataoutput.write( '' ) dataoutput.write( 'Distribution of events by lunar phases: \n' ) for elem in fasi: dataoutput.write( elem + ', ' + str(fasi[elem]) + '\n') dataoutput.write( '' ) #dof = len(fasi)-1 #distance =chi2UniformDistance(fasi) #prob = chi2Probability( dof, distance) #dataoutput.write( "Chi Square: dof: %d, distance: %.4f, P-value: %.4f \n\n" % (dof, distance, prob) ) dataoutput.write("\n\n") dataoutput.write(":::: Perigees ::::\n") perigees.pop(0) somma = 0 for perigeo in perigees: somma = somma + perigeo[1] dataoutput.write(str(perigeo[0]) + ',' + str(perigeo[1]) + '\n') #chi, pvalue = chisquare(perigees) #dataoutput.write( "Chi Square: %.3f, P-value: %.4f \n" % (chi, pvalue) ) media = somma / len(perigees) dataoutput.write( "Mean No. of Events in Perigees: %.2f \n\n" % media ) dataoutput.write(":::: Apogees ::::\n") somma = 0 apogees.pop(0) for apogeo in apogees: somma = somma + apogeo[1] dataoutput.write(str(apogeo[0]) + ',' + str(apogeo[1]) + '\n') dataoutput.write("\n\n") #chi, pvalue = chisquare(apogees) #dataoutput.write( "Chi Square: %3f, P-value: %.4f \n" % (chi, pvalue) ) media = somma / len(apogees) dataoutput.write("Mean No. of Events in Apogees: %.2f \n\n" % media ) dataoutput.write(":::: phases ::::\n") max = 0 for item in sinodicofasi: dataoutput.write(item + ", ") sinodicofasi[item][0] = len(sinodicofasi[item]) if sinodicofasi[item][0] >= max: max = sinodicofasi[item][0] dataoutput.write("\n") for j in range(1, (max)): for phase in sinodicofasi: if j < sinodicofasi[phase][0]: dataoutput.write(str(sinodicofasi[phase][j]) + ",") else: dataoutput.write(",") dataoutput.write("\n") dataoutput.write("\n") dataoutput.close() print 'File of output completed, see file %s' % (path.strip('./')) print '' print 'Now writing the matrix events / sinodic month' path = './matricesinodica' + timestamp + '.txt' dataoutput = open(path,'a') dataoutput.write(":::: vertical ::::\n") max = 0 for elem in sinodico: elem[0] = len(elem) if elem[0] >= max: max = elem[0] for j in range(1, (max)): for i in range(0, len(sinodico)): if j < sinodico[i][0]: dataoutput.write(str(sinodico[i][j]) + ",") else: dataoutput.write(",") dataoutput.write("\n") dataoutput.write("\n") dataoutput.write(":::: horizontal ::::\n") for elem in sinodico: elem.pop(0) #elimino il primo elemento, che era solo un placeholder for j in range(0, len(elem)): dataoutput.write(str(elem[j]) + ", ") dataoutput.write("\n") dataoutput.close() print 'Matrix events / sinodic month completed, see file %s' % (path.strip('./')) print '' print 'Now writing the temporal series file from %s to %s' % (datesince, dateuntil) path = './serietemporale' + timestamp + '.txt' fileoutput = open(path,'a') cicloanno = datesince.strftime("%Y") ciclomese = datesince.strftime("%m") ciclogiorno = datesince.strftime("%d") data = datesince prog = 0 while data <= dateuntil: datascritta = str(cicloanno) + "-" + str(ciclomese) + "-" + str(ciclogiorno) numeronascite = contenutofile.count(datascritta) #serietemporale.append(numeronascite) #serietemporale[indice] = numeronascite elencoperigei = str(perigees) elencoapogei = str(perigees) stringaadd = '' if elencoperigei.count(datascritta) > 0: stringaadd = ',P,' elif elencoapogei.count(datascritta) > 0: stringaadd = ',A,' fileoutput.write(datascritta + ',' + str(numeronascite) + stringaadd + "\n") data = data + RelativeDateTime(days=+1) cicloanno = data.strftime("%Y") ciclomese = data.strftime("%m") ciclogiorno = data.strftime("%d") prog +=1 if (prog % 30) == 0: print 'Now written value for date %s-%s-%s' % (cicloanno, ciclomese, ciclogiorno) fileoutput.close() print 'File of the Temporal Series completed, see file %s' % (path.strip('./')) print '' print 'Now writing a list of the plenilunes and novilunes found in the considered interval' path = './PleniNoviList' + timestamp + '.txt' fileoutput = open(path,'a') fileoutput.write('Pleniluni \n' + elencopleni + '\n' + 'Noviluni\n' + elenconovi) fileoutput.close() print 'File of the plenilune, novilune list completed, see file %s' % (path.strip('./')) print ' ******* Exit from script *******' if __name__ == "__main__": main()