mirror of
https://asciireactor.com/otho/psdlag-agn.git
synced 2024-11-22 10:35:08 +00:00
231 lines
6.8 KiB
C++
231 lines
6.8 KiB
C++
/*
|
|
* lag.cpp
|
|
*
|
|
* Created on: Jun 1, 2013
|
|
* Author: azoghbi
|
|
*/
|
|
|
|
#include "inc/lag_rms.hpp"
|
|
|
|
lag_rms::lag_rms( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
|
|
|
|
// ----------- initial parameters ------------ //
|
|
n1 = lc1.len;
|
|
n = n1 + lc2.len;
|
|
dt = lc1.dt;
|
|
// ------------------------------------------ //
|
|
|
|
|
|
// ----------- light curve setup ------------ //
|
|
setlc();
|
|
mean1 = lc1.demean(); mean2 = lc2.demean();
|
|
mean1sq = mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
|
|
int i;
|
|
for( i=0 ; i<n1 ; i++ ){
|
|
lc[i] = lc1.lc[i];
|
|
lce[i] = lc1.lce[i]*lc1.lce[i];
|
|
t[i] = lc1.t[i];
|
|
}
|
|
for( i=n1 ; i<n ; i++ ){
|
|
lc[i] = lc2.lc[i-n1];
|
|
lce[i] = lc2.lce[i-n1]*lc2.lce[i-n1];
|
|
t[i] = lc2.t[i-n1];
|
|
}
|
|
// ------------------------------------------ //
|
|
|
|
// ----------- parent initilizer ------------ //
|
|
init( fqL , 2 );
|
|
// ------------------------------------------ //
|
|
|
|
// --------- constants and indices ---------- //
|
|
p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
|
|
for( i=0 ; i<nfq ; i++ ){
|
|
p1[i] = pars[i] * mean1sq;
|
|
p2[i] = pars[i+nfq] * mean2sq;
|
|
icx[i] = i;
|
|
iphi[i] = i+nfq;
|
|
}
|
|
// ------------------------------------------ //
|
|
|
|
}
|
|
|
|
lag_rms::~lag_rms() {}
|
|
|
|
void lag_rms::_C( vec p ){
|
|
int i,j,k; double d;
|
|
for( i=0 ; i<nfq ; i++ ){ p[icx[i]] *= mean12;}
|
|
for( i=0 ; i<n1 ; i++ ){
|
|
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
|
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); } C(i,j) = d;}
|
|
}
|
|
for( i=n1 ; i<n ; i++ ){
|
|
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); } d += lce[i]; C(i,i) = d;
|
|
for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); } C(i,j) = d;}
|
|
for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
|
|
d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); } C(i,j) = d;}
|
|
}
|
|
}
|
|
|
|
void lag_rms::_dC( vec p , int k ){
|
|
if( k<nfq ){ _dCx( p , k );
|
|
}else{ _dPhi( p , k-nfq );}
|
|
}
|
|
|
|
void lag_rms::_dCx( vec p , int k ){
|
|
int i,j; double phi=p[iphi[k]];
|
|
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
|
for( i=n1 ; i<n ; i++ ){
|
|
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
|
for( j=0 ; j<n1 ; j++ ){
|
|
dC(i,j) = ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
|
|
}
|
|
}
|
|
|
|
void lag_rms::_dPhi( vec p , int k ){
|
|
int i,j; double cx=p[k]*mean12,phi=p[iphi[k]];
|
|
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
|
for( i=n1 ; i<n ; i++ ){
|
|
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
|
for( j=0 ; j<n1 ; j++ ){
|
|
dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
|
|
}
|
|
}
|
|
|
|
|
|
void lag_rms::step_pars( int n , vec& dpar , vec& pars ){
|
|
for( int i=0 ; i<npar ; i++ ){
|
|
if( dpar[i]>3000 ){dpar[i] = 3000;} if( dpar[i]<-3000 ){dpar[i] = -3000;}
|
|
pars[i] += dpar[i]/((n<10)?10:1);
|
|
}
|
|
}
|
|
|
|
|
|
void lag_rms::print_pars( vec& pars , vec& errs ){
|
|
for( int i=0 ; i<nfq ; i++ ){
|
|
if( pars[i]<0 ){ pars[i]*=-1; pars[i+nfq] += M_PI;}
|
|
while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
|
|
while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
|
|
}
|
|
mod::print_pars( pars , errs );
|
|
}
|
|
|
|
|
|
// ****************************************** //
|
|
|
|
|
|
lag10_rms::lag10_rms( lcurve lc1 , lcurve lc2 , vec fqL , vec pars ) {
|
|
|
|
// ----------- initial parameters ------------ //
|
|
n1 = lc1.len;
|
|
n = n1 + lc2.len;
|
|
dt = lc1.dt;
|
|
// ------------------------------------------ //
|
|
|
|
|
|
// ----------- light curve setup ------------ //
|
|
setlc();
|
|
mean1 = lc1.demean(); mean2 = lc2.demean();
|
|
mean1sq = mean1*mean1; mean2sq = mean2*mean2; mean12 = mean1*mean2;
|
|
int i;
|
|
for( i=0 ; i<n1 ; i++ ){
|
|
lc[i] = lc1.lc[i];
|
|
lce[i] = lc1.lce[i]*lc1.lce[i];
|
|
t[i] = lc1.t[i];
|
|
}
|
|
for( i=n1 ; i<n ; i++ ){
|
|
lc[i] = lc2.lc[i-n1];
|
|
lce[i] = lc2.lce[i-n1]*lc2.lce[i-n1];
|
|
t[i] = lc2.t[i-n1];
|
|
}
|
|
// ------------------------------------------ //
|
|
|
|
// ----------- parent initilizer ------------ //
|
|
init( fqL , 2 );
|
|
// ------------------------------------------ //
|
|
|
|
// --------- constants and indices ---------- //
|
|
p1.setlength(nfq); p2.setlength(nfq); icx.setlength(nfq); iphi.setlength(nfq);
|
|
for( i=0 ; i<nfq ; i++ ){
|
|
p1[i] = pow(10,pars[i]) * mean1sq;
|
|
p2[i] = pow(10,pars[i+nfq]) * mean2sq;
|
|
icx[i] = i;
|
|
iphi[i] = i+nfq;
|
|
}
|
|
// ------------------------------------------ //
|
|
|
|
}
|
|
|
|
lag10_rms::~lag10_rms() {}
|
|
|
|
void lag10_rms::_C( vec p ){
|
|
int i,j,k; double d; for(i=0;i<nfq;i++){p[icx[i]] = pow(10,p[i]) * mean12;}
|
|
for( i=0 ; i<n1 ; i++ ){
|
|
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,i); }
|
|
d += f1*p1[0]*Cfq2[0](i,i) + f2*p1[nfq-1]*Cfq2[2](i,i);
|
|
d += lce[i]; C(i,i) = d;
|
|
for( j=0 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p1[k]*Cfq[k](i,j); }
|
|
d += f1*p1[0]*Cfq2[0](i,j) + f2*p1[nfq-1]*Cfq2[2](i,j);
|
|
C(i,j) = d;}
|
|
}
|
|
for( i=n1 ; i<n ; i++ ){
|
|
d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,i); }
|
|
d += f1*p2[0]*Cfq2[0](i,i) + f2*p2[nfq-1]*Cfq2[2](i,i);
|
|
d += lce[i]; C(i,i) = d;
|
|
for( j=n1 ; j<i ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){ d += p2[k]*Cfq[k](i,j); }
|
|
d += f1*p2[0]*Cfq2[0](i,j) + f2*p2[nfq-1]*Cfq2[2](i,j);
|
|
C(i,j) = d;}
|
|
for( j=0 ; j<n1 ; j++ ){ d = 0; for( k=0 ; k<nfq ; k++ ){
|
|
d += p[icx[k]] * ( Cfq[k](i,j)*cos(p[iphi[k]]) - Sfq[k](i,j)*sin(p[iphi[k]])); }
|
|
d += f1*p[icx[0]] * ( Cfq2[0](i,j)*cos(p[iphi[0]]) - Sfq2[0](i,j)*sin(p[iphi[0]]) );
|
|
d += f2*p[icx[nfq-1]] * ( Cfq2[2](i,j)*cos(p[iphi[nfq-1]]) - Sfq2[2](i,j)*sin(p[iphi[nfq-1]]) );
|
|
C(i,j) = d;}
|
|
}
|
|
}
|
|
|
|
void lag10_rms::_dC( vec p , int k ){
|
|
if( k<nfq ){ _dCx( p , k );
|
|
}else{ _dPhi( p , k-nfq );}
|
|
}
|
|
|
|
void lag10_rms::_dCx( vec p , int k ){
|
|
int i,j; double phi=p[iphi[k]],cx = log(10)*pow(10,p[k])*mean12;
|
|
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
|
for( i=n1 ; i<n ; i++ ){
|
|
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
|
for( j=0 ; j<n1 ; j++ ){
|
|
dC(i,j) = cx * ( Cfq[k](i,j)*cos(phi) - Sfq[k](i,j)*sin(phi) ); dC(j,i) = dC(i,j);}
|
|
}
|
|
}
|
|
|
|
void lag10_rms::_dPhi( vec p , int k ){
|
|
int i,j; double cx=pow(10,p[k])*mean12,phi=p[iphi[k]];
|
|
for( i=0 ; i<n1 ; i++ ){ for( j=0 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;} }
|
|
for( i=n1 ; i<n ; i++ ){
|
|
for( j=n1 ; j<=i ; j++ ){ dC(i,j) = 0; dC(j,i) = 0;}
|
|
for( j=0 ; j<n1 ; j++ ){
|
|
dC(i,j) = cx * ( -Cfq[k](i,j)*sin(phi) - Sfq[k](i,j)*cos(phi) ); dC(j,i) = dC(i,j);}
|
|
}
|
|
}
|
|
|
|
|
|
void lag10_rms::step_pars( int n , vec& dpar , vec& pars ){
|
|
for( int i=0 ; i<npar ; i++ ){
|
|
if( dpar[i]>3 ){dpar[i] = 3;} if( dpar[i]<-3 ){dpar[i] = -3;}
|
|
//pars[i] += dpar[i];
|
|
pars[i] += dpar[i]/((n<5)?10:1);
|
|
}
|
|
}
|
|
|
|
|
|
void lag10_rms::print_pars( vec& pars , vec& errs ){
|
|
for( int i=0 ; i<nfq ; i++ ){
|
|
while( pars[i+nfq] > M_PI ){ pars[i+nfq] -= 2*M_PI; }
|
|
while( pars[i+nfq] <-M_PI ){ pars[i+nfq] += 2*M_PI; }
|
|
}
|
|
mod::print_pars( pars , errs );
|
|
}
|
|
|
|
void lag10_rms::what_pars( int& ip1 , int& ip2 ){
|
|
ip1 = nfq; ip2 = npar;
|
|
}
|