Hi Doug

I use IntegrationTools:-Expand because I was told to follow that way here in the forum! This was becaus of the symbols in the integrand.

I can copy a sheet of the code, till the scalar product functions. Thank for any help

restart;

time_start:= time():

Digits:=50:

S_h:=proc(t)

local f_0,T;

global S_0:

#S_0:=1.0*10^(-49):

f_0 :=215.0:

T:=t/f_0:

simplify(S_0*(T^(-4.14) -5*T^(-2)+ 111*(1-T^2+T^4/2)/(1+T^2/2)))

end proc:

assume(eta,'positive');

assume(phi,'positive');

assume(t,'positive');

assume(Enn,'positive');

assume(M,'positive');

assume(f,'positive');

assume(Emm,'positive');

assume(S_0,'positive');

assume(rho,'positive');

with(LinearAlgebra):

h1_digits:=50;

h2_digits:=50;

h3_digits:=50;

parameters:= [t,phi,Emm,eta]:

param_dim:= nops(parameters):

f_low:=20.0:

M_solar:=4.927*10^(-6):

eta_eff:= .25:

PNorder:=2:

M:=0.0;

`index/commutingindexes`:=proc(idx,M,val)

local idx1:

idx1:=op(sort(idx));

if nargs = 2 then

# retrieving

M[idx1];

else

# storing

M[idx1]:=op(val);

end if:

end proc:

for counter from 0 by 1 to 19 do ##################### FOR CYCLE

M_tot:=(2.8 + 219.0*counter/20)*M_solar:

appendto("./M_tot_seq.txt"):

lprint(M_tot);

appendto("./times.txt");

printf("Counter value is %d \n", counter);

appendto("./garbage.txt"):

f_lso:= (6^(3/2)*Pi*M_tot)^(-1):

v_lso:=(Pi*M_tot*f_lso)^(1/3):

f_high:=f_lso:

variables:=[f]:

variab_dim:=nops(variables):

p:= proc()

option hfloat;

local Psi,i,psi,alpha_par,v,theta,euler_gamma,lambda;

global Enn,PNorder,param_dim,f,eta,M,M_tot,f_lso,v_lso,Emm,eta_eff,parameters,variables,S_0,rho,t,phi;

Digits:=50;

M:=Emm*eta^(-3/5);

v:=(Pi*M*f)^(1/3);

theta:=-1.28;

lambda:=-0.6451;

euler_gamma:=0.577215664901532860606512090082:

alpha_par:=[

1

,0

,20/9*(743/336 + 11/4*eta)

,-16*Pi

,10*(3058673/1016064+5429/1008*eta+617/144*eta^2)

, Pi*(38645/756+38645/252*ln(v/v_lso) -65/9*eta*(1+3*ln(v/v_lso)))

, (11583231236531/4694215680-640/3*Pi^2-6848/21*euler_gamma)+eta*(-15335597827/3048192+2255/12*Pi^2-1760/3*theta+12320/9*lambda)+76055/1728*eta^2-127825/1296*eta^3-6848/21*ln(4*v)

,Pi*(77096675/254016+378515/1512*eta-74045/756*eta^2)];

Psi:= (2*Pi*f*t-phi-Pi/4+ 3/(128*eta*v^5)*add(alpha_par[i+1]*v^i, i=0..PNorder));

Enn*f^(-7/6)*exp(I*Psi):

end proc:

H1:=proc()

global Enn,PNorder,param_dim,f,eta,M,M_tot,f_lso,v_lso,Emm,eta_eff,parameters,variables,S_0,rho,t,phi,h1, counter_h1,h1_digits;

Digits:=h1_digits;

h1:=Array(1..param_dim, i->

'eval'(

(diff(p(), parameters[i]))

,[Emm=(M_tot*(eta_eff)^(3/5)),eta=eta_eff] )

):

end proc;

time_h1:=time(H1());

appendto("./time-h.txt"):

printf("h1 has been calculated in %a seconds with %d digits \n", time_h1, h1_digits);

appendto("./garbage.txt"):

prova:=simplify(h1):

save prova,"./h1":

scalar_pr:= proc(A,B)

option hfloat;

global variables, parameters,f_low, f_high,S_h,M_tot,eta_eff,S_0,f_lso,f,Enn,rho;

local num;

description "Calculate the scalar product between two tensors A and B in the frequency space":

Digits:=15;

num:=

Re(

simplify(

A*conjugate(B)+conjugate(A)*B

)

)

;

evalf(

IntegrationTools:-Expand(

2.0*Int(num/S_h(f), f=f_low..f_high)

)

);

appendto("./time_snr.txt");

lprint(time_snr_end);

appendto("./garbage.txt"):

Low_fisher:=proc()

option hfloat;

description " Calculate the fisher information matrix called low_fisher(a,b)";

global Enn,low_fisher,S_0, f_low,f_high,param_dim,M_tot,eta_eff,parameters,variables,rho,t,phi,M,eta,Emm:

local gg,a_1,b_1:

gg:=(a_1,b_1)->

scalar_pr(h1[a_1],h1[b_1])

;

low_fisher:= Matrix(param_dim,param_dim, shape=symmetric,gg);

end proc:

Low_fisher():

#save low_fisher, "./low_fisher.txt":

Up_fisher:=proc()

option hfloat;

description " Calculate the upper index fisher information matrix called up_fisher(a,b)";

global Enn, up_fisher,M_tot,eta_eff,S_0,param_dim,variables,parameters,rho,t,phi,Emm,eta,M:

local a_3,b_3,gg;

gg:=(a_3,b_3)-> (MatrixInverse(low_fisher))[a_3,b_3];

up_fisher:= Matrix(param_dim,param_dim,shape=symmetric,gg);

end proc:

Up_fisher():

save up_fisher,"./up_fisher":

end proc:

time_snr_in:=time():

snr:=scalar_pr(p(),p()):

time_snr_end:=time()-time_snr_in:

save snr, "./snr_squared.txt":

appendto("./snr-seq.txt"):

snr;

end do;