{-

Copyright 2015 Ken Takusagawa

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.

-}
{-# LANGUAGE ScopedTypeVariables #-}
module Main where {
import qualified Data.List as List;
import qualified Data.Ord as Ord;
import Text.Printf(printf);
a_freq :: Double;
a_freq = 440;
-- also consider 44100 and 48000 and 1

fundamental_frequency :: Double;
fundamental_name :: String;
fundamental_frequency = planck_frequency;
fundamental_name = "Pf";

planck_frequency :: Double;
planck_frequency = let {
plancktime :: Double;
-- plancktime = 5.391232e-44;  --from units
plancktime = 5.3910604e-44; -- via definition and G
} in 1/ plancktime;

hydrogen_frequency :: Double;
hydrogen_frequency = 1420.40575177e6;

cesium_frequency :: Double;
cesium_frequency = 9192631770;

equal_temperament :: Double;
equal_temperament = 2 ** (1/12);

year :: Double;
-- gregorian year in seconds
year = 60 * 60 * 24 * 365.2425;
-- piano range
lo :: Double;
--lo = a_freq*equal_temperament**(negate 48);
lo = 1/ (150 * year);
hi :: Double;
--hi = a_freq*equal_temperament**(39);
hi = 1/year;

freq_to_note :: Double -> Double;
freq_to_note f = log(f/a_freq)/log equal_temperament;

cents_calc :: Double -> Double;
cents_calc f = freq_to_note f - (fromIntegral ((round $ freq_to_note f)::Integer));

do_base :: Integer -> [(Double,Metadata)];
do_base base = do {
 let { logbase :: Double -> Double; logbase x = log(x)/ (log $ fromIntegral base); };
 i :: Integer <- enumFromTo
  (ceiling  $ logbase (fundamental_frequency / hi) )
  (floor $ logbase (fundamental_frequency / lo) );
  let {
   f = getfreq base i;
   n = round $ freq_to_note $ getfreq base i;
   c = cents_calc f;
   (octave,remainder) = divMod (n+9) 12; -- 9 is the semitone position of A in the C scale
  };
  return (100*c, (base, i, n, (List.genericIndex note_names remainder, octave+4),f, a_adjusted c ));
};

a_adjusted :: Double -> Double;
a_adjusted dx = equal_temperament**dx *a_freq;

type Metadata = (Integer,Integer,Integer,(String,Integer),Double,Double);
all_bases :: Integer -> [(Double,Metadata)];
all_bases n = concat $ map do_base $ [2..n];

sorted_bases :: Integer -> [(Double,Metadata)];
sorted_bases = List.sortBy (Ord.comparing (abs . fst)) . all_bases;

zipMap :: (a -> b) -> [a] -> [(a,b)];
zipMap f l = zip l $ map f l;

getfreq :: Integer -> Integer -> Double;
getfreq base i = fundamental_frequency / (fromIntegral $ base^i);

examine :: Integer -> IO();
examine = mapM_ print . take 20 . sorted_bases;

base_of :: Metadata -> Integer;
base_of (i,_,_,_,_,_) = i;

-- best approximation with a base less than or equal to starting_base
chase_down :: Integer -> [(Double,Metadata)];
chase_down 0 = [];
chase_down starting_base = case sorted_bases starting_base of {
[] -> [];
(h:_) -> h:(chase_down $ pred $ base_of $ snd h)
};

note_names :: [String];
note_names = ["C","C#","D","D#","E","F","F#","G","G#","A","A#","B"];

output_line :: (Double, Metadata) -> String;
output_line (cents,(base,expo,_notenum,(note,octave),freq,avalue)) =
 note ++ show octave ++ " = " ++ fundamental_name ++ " / " ++ show base ++ "^" ++ show expo ++ " = " ++ (printf "%.4f" freq) ++ " Hz, or A = " ++ (printf "%.4f" avalue) ++ " Hz, offset = " ++ (printf "%.5f" cents) ++ " cents";

main :: IO();
main = do{
--  print a_freq;
--  mapM_ (putStrLn . output_line) $ chase_down 1000000;
    go_year_tuples;
};
-- examine 100000;
-- 40000 gets everything until 3899672
{-
selected entries
planckfrequency
1.8548636007502552e43

C3 (-0.00005282066055656287,(37360,9,-21,(-2,3),130.81277865914763,439.9999865754198))
37360: 2 2 2 2 5 467

A#2 (0.00010499820106701918,(5567,11,-23,(-2,1),116.54094744765361,440.0000266857101))
5567: 19 293

F#4 (-0.0194702213589526,(5012,11,-3,(-1,9),369.9902616146372,439.9950515938448))
5012: 2 2 7 179

A1 (0.038827264778262816,(394,16,-36,(-3,0),55.00123352674986,440.00986821399965))
394: 2 197

F#3 (4.571492368973473,(10,41,-15,(-2,9),185.48636007502552,441.16339827435434))

A7 (8.616156778325035,(7,47,36,(3,0),3537.5623099173445,442.19528873966726))

digression, two new units of time
? t*10^43
%55 = 0.5391232000000000000000000000
? t*2^143
%56 = 0.6011424556908232423029394675

-- and now using the new value of plancktime

D#5 (5.331509575157156e-4,(31416,9,6,(0,6),622.2541590735664,440.00013550245126))

C#5 (1.176172096766237e-2,(1310,13,4,(0,4),554.3690282334554,440.00298929818774))

A#3 (-9.058358170861425e-2,(360,16,-11,(-1,1),233.0696855048837,439.9769784257213))

A#1 (0.3241916026510694,(77,22,-35,(-3,1),58.28138296379532,440.0824022967275))

D#4 (0.43060043597407116,(50,24,-6,(-1,6),311.2043782703677,440.10945241984103))

E1 (-0.667694537656871,(40,26,-41,(-4,7),41.187556519498706,439.8303355048625))

B5 (3.9505950053431604,(22,30,14,(1,2),990.0232133910473,441.0052058696316))
(this is very close to A441, not exceeded until)
(2.4035428722868346e-3,(6022,11,-38,(-4,10),49.11086002004969,441.0006122587198))
(2.043656590977605e-2,(22,30,14,(1,2),990.0232133910473,441.0052058696316))

F#3 (4.626597489818707,(10,41,-15,(-2,9),185.49226419351564,441.1774407137868))

A7 (8.671261899171157,(7,47,36,(3,0),3537.674912197294,442.20936402466117))

-- A442
(-1.9391513063737875e-3,(30214,9,12,(1,0),883.999009834173,441.99950491708626))
(2.416789365566885e-3,(11744,10,-3,(-1,9),371.67673439963113,442.00061702901803))
(-7.6909594618257415e-3,(217,17,36,(3,0),3535.98429145253,441.99803643156565))
(9.394296122167134e-3,(151,19,-31,(-3,5),73.75030213736711,442.0023984568431))
(-0.18071576983516024,(62,23,-24,(-2,0),110.48846600182573,441.95386400730337))
(0.29903146805274616,(38,26,-18,(-2,6),156.2975931797607,442.0763520821806))
(0.6924432846861173,(37,26,-6,(-1,6),312.6662296793513,442.1768225085999))
(0.8198468590443042,(7,47,36,(3,0),3537.674912197294,442.20936402466106))

-- A443
(4.557859890041982e-3,(25564,9,38,(3,2),3978.01597220922,443.0011662978731))
(-7.2240594832351235e-3,(5988,11,-37,(-4,11),52.266821966324116,442.99815146214587))
(4.894761045335372e-2,(391,16,-34,(-3,2),62.158093311505304,443.012525226142))
(7.422264178345017e-2,(226,17,24,(2,0),1772.0759720158671,443.0189930039663))
(-0.8030754796628514,(163,18,32,(2,8),2811.5701456041993,442.79455121283985))
(1.5446171165041989,(60,23,-11,(-1,1),234.88054373338247,443.3954235572978))
(-1.9050727713402438,(35,26,19,(1,7),1326.0400749807463,442.51278507772986))
(-2.0710005993798575,(34,27,-29,(-3,7),82.8695617825808,442.47037505997974))
(2.6949987960477984,(18,33,-32,(-3,4),69.87682022722807,443.6901517475791))
(-3.0925481372591435,(7,47,36,(3,0),3537.674912197294,442.2093640246612))

Inspired by Scientific Pitch C256, A = 430.5389646099018460319362438

-}

year_tuples :: [(Integer,Integer)];
year_tuples = do {(_,(x,y,_,_,_,_))<-all_bases 100; return (x,y)};

display_year_tuple :: (Integer,Integer) -> String;
display_year_tuple (x,y) = unwords [show $ (fromIntegral (x^y)/planck_frequency/year),show x, show y];

go_year_tuples :: IO ();
go_year_tuples = mapM_ (putStrLn . display_year_tuple) year_tuples;

} -- end