# génération des courbes de réponse en fréquence des corrections Phono
# copyright Pierre Ragnarsson
# Programme et fichiers générés utilisables seulement à des fins personnelles et non commerciales
# v1 du 16 Octobre 2023

import math
import cmath
import numpy as np

f1=20                        # fréquence inférieure
f2=20000                     # fréquence supérieure
N=1000                       # nombre de points
offset1=20.0                   # décalage en dB pour les courbes permettant de prendre en compte le gain
offset2=-20.0                  # décalage en dB pour les courbes inverse permettant de prendre en compte le gain
t1=0.00318                   # t1=3180us	
t2=0.000318                  # t2=318us	
t3=0.000075                  # t3=75us
t4=0.00000296736             # t4=2.967us (337kHz)	
t5=0.00795                   # t5=7950us


file_riaa=open("riaa.frd","w")               # nom du fichier .frd contenant la réponse en fréquence RIAA
file_inv_riaa=open("inv_riaa.frd","w")       # nom du fichier .frd contenant la réponse en fréquence inverse RIAA
file_iec=open("riaa_iec.frd","w")            # nom du fichier .frd contenant la réponse en fréquence RIAA avec amendement IEC
file_inv_iec=open("inv_riaa_iec.frd","w")    # nom du fichier .frd contenant la réponse en fréquence inverse RIAA avec amendement IEC
       
for i in range(0,N):
    f=f1*(f2/f1)**(i/(N-1))
    
    # calcul des fonctions de transfert
    riaa1=2*np.pi*f*t1
    riaa2=2*np.pi*f*t2
    riaa3=2*np.pi*f*t3
    riaa4=2*np.pi*f*t4
    iec_amd=2*np.pi*f*t5
    inv_riaa=((1+1j*riaa1)*(1+1j*riaa3))/((1+1j*riaa2)*(1+1j*riaa4))
    inv_iec=inv_riaa*(1+1j*iec_amd)/(1j*iec_amd)
    riaa=((1+1j*riaa2)*(1+1j*riaa4))/((1+1j*riaa1)*(1+1j*riaa3))
    iec=riaa*(1j*iec_amd)/(1+1j*iec_amd)
    
    #inv_riaa=1.0/riaa
    #iec=1.0/inv_iec
    
    file_riaa.write(str(f))
    file_riaa.write("\t")
    file_riaa.write(str(20.0*math.log10(abs(riaa))+offset1))
    file_riaa.write("\t")
    file_riaa.write(str((180.0/np.pi)*cmath.phase(riaa)))
    file_riaa.write("\n")

    file_inv_riaa.write(str(f))
    file_inv_riaa.write("\t")
    file_inv_riaa.write(str(20.0*math.log10(abs(inv_riaa))+offset2))
    file_inv_riaa.write("\t")
    file_inv_riaa.write(str((180.0/np.pi)*cmath.phase(inv_riaa)))
    file_inv_riaa.write("\n")

    file_iec.write(str(f))
    file_iec.write("\t")
    file_iec.write(str(20.0*math.log10(abs(iec))+offset1))
    file_iec.write("\t")
    file_iec.write(str((180.0/np.pi)*cmath.phase(iec)))
    file_iec.write("\n")
    
    file_inv_iec.write(str(f))
    file_inv_iec.write("\t")
    file_inv_iec.write(str(20.0*math.log10(abs(inv_iec))+offset2))
    file_inv_iec.write("\t")
    file_inv_iec.write(str((180.0/np.pi)*cmath.phase(inv_iec)))
    file_inv_iec.write("\n")
    
file_riaa.close()
file_inv_riaa.close()
file_iec.close()
file_inv_iec.close()