pro p_phs, ps, in, res1, res2, par1, par2, ndeg = ndeg
;
; This is to fit the MS phase-intensity plot for different periods in order
l to derive the ratio between the old and new filters.  
;

ss1 = where(ps gt 0 and in gt 10)
ss2 = where(ps lt 0 and in gt 10)
help, ss1, ss2
ps1 = abs(ps(ss1))
ps2 = abs(ps(ss2))
in1 = in(ss1)
in2 = in(ss2)
ord1 = sort(ps1)
ord2 = sort(ps2)
ps1 = ps1(ord1)
ps2 = ps2(ord2)
in1 = in1(ord1)
in2 = in2(ord2)
if n_elements(ndeg) eq 0 then ndeg = 3
par1 = poly_fit(abs(ps1), in1, ndeg, yfit = yfit1)
par2 = poly_fit(abs(ps2), in2, ndeg, yfit = yfit2)

res1 = fltarr(181)
res2 = fltarr(181)
ph = findgen(181)
for j = 0, ndeg do res1 = res1 + (abs(ph))^(j(0))*par1(j)
for j = 0, ndeg do res2 = res2 + (abs(ph))^(j(0))*par2(j)

key = ''
!p.multi = [0, 0, 2, 0, 0]
plot, ps1, in1, psym = 1, color = 200
oplot, res1, psym = 1
plot, ps2, in2, psym = 1, color = 200
oplot, res2, psym = 1
read, key
end

	pro check_2period, ps0, ps1, in0, in1, cin0, cin1, ein0, ein1, ndeg = ndeg

;
;
;

;restore, 'new_old_period_pahse_int.dat'

!p.charsize = 1.5
p_phs, ps0, cin0, cres01, cres02, ndeg = ndeg
p_phs, ps1, cin1, cres11, cres12, ndeg = ndeg
 
p_phs, ps0, in0, mres01, mres02, ndeg = ndeg
p_phs, ps1, in1, mres11, mres12, ndeg = ndeg
 
p_phs, ps0, ein0, eres01, eres02, ndeg = ndeg
p_phs, ps1, ein1, eres11, eres12, ndeg = ndeg

set_plot, 'ps'
device, filename = 'check_2period.ps', xs = 18, ys = 22, yoff = 2, /color
!p.multi = [0, 0, 3, 0, 0]
plot, ps0, alog(cin0), /ynoz, psym = 1, tit = 'CS', xr = [-140, 140], xst = 1, yr = [31, 36], yst = 1
oplot, ps1, alog(cin1), color = 200, psym = 1
plot, ps0, alog(in0), /ynoz, psym = 1, tit = 'MS', xr = [-140, 140], xst = 1, yr = [22, 24.2], yst = 1
oplot, ps1, alog(in1), color = 200, psym = 1
plot, ps0, alog(ein0>1.), /ynoz, psym = 1, tit = 'ES', xr = [-140, 140], xst = 1, yr = [20, 23], yst = 1
oplot, ps1, alog(ein1>1.), color = 200, psym = 1

!p.multi = [0, 0, 3, 0, 0]

loadct, 38
plot, [-(indgen(181)), indgen(181)], [alog(cres02>1.), alog(cres01>1.)], psym = 1, xr = [-140, 140], xst = 1, yr = [31, 36], yst = 1
oplot, (ps1), alog(cin1>1.), color = 200, psym = 1
 
plot, [-(indgen(181)), indgen(181)], [alog(mres02>1.), alog(mres01>1.)], psym = 1, xr = [-140, 140], xst = 1, yr = [22, 24.2], yst = 1
oplot, (ps1), alog(in1>1.), color = 200, psym = 1
 
plot, [-(indgen(181)), indgen(181)], [alog(eres02>1.), alog(eres01>1.)], psym = 1, xr = [-140, 140], xst = 1, yr = [20, 23], yst = 1
oplot, (ps1), alog(ein1>1.), color = 200, psym = 1
 

s1 = where(ps1 gt 0 and in1 gt 0)
s2 = where(ps1 lt 0 and in1 gt 0)
ord1 = sort(abs(ps1(s1)))
ord2 = sort(abs(ps1(s2)))
ps11 = (abs(ps1(s1)))(ord1)
ps12 = (abs(ps1(s2)))(ord2)
in11 = (in1(s1))(ord1)
in12 = (in1(s2))(ord2)
cin11 = (cin1(s1))(ord1)
cin12 = (cin1(s2))(ord2)
ein11 = (ein1(s1))(ord1)
ein12 = (ein1(s2))(ord2)
 
cs1 = spline(indgen(181), cres01, ps11)
cs2 = spline(indgen(181), cres02, ps12)
ms1 = spline(indgen(181), mres01, ps11)
ms2 = spline(indgen(181), mres02, ps12)
es1 = spline(indgen(181), eres01, ps11)
es2 = spline(indgen(181), eres02, ps12)

help, cin11, cin12, in11, in12, ein11, ein12
help, cs1, cs2, ms1, ms2, es1, es2

!p.multi = [0, 0, 3, 0, 0]
plot, [-ps12, ps11], [cs2, cs1], psym = 1
oplot, ps11, cin11, psym = 1, color = 200  
oplot, -ps12, cin12, psym = 1, color = 200  
plot, [-ps12, ps11], [ms2, ms1], psym = 1
oplot, ps11, in11, psym = 1, color = 200  
oplot, -ps12, in12, psym = 1, color = 200  
plot, [-ps12, ps11], [es2, es1], psym = 1
oplot, ps11, ein11, psym = 1, color = 200  
oplot, -ps12, ein12, psym = 1, color = 200  

device, /close
set_plot,'x'
print, sqrt(total(((cs1 - cin11)/cs1)^2))/(n_elements(cs1)), mean((cs1-cin11)/cs1)
print, sqrt(total(((cs2 - cin12)/cs2)^2))/(n_elements(cs2)), mean((cs2-cin12)/cs2)
print, sqrt(total(((ms1 - in11)/ms1)^2))/(n_elements(ms1)), mean((ms1-in11)/ms1)
print, sqrt(total(((ms2 - in12)/ms2)^2))/(n_elements(ms2)),mean((ms2-in12)/ms2)
print, sqrt(total(((es1 - ein11)/es1)^2))/(n_elements(es1)),mean((es1-ein11)/es1)
print, sqrt(total(((es2 - ein12)/es2)^2))/(n_elements(es2)),mean((es2-ein12)/es2)

print, ' '
print, mean(cres01(60:120))/mean(cres11(60:120))
print, mean(cres02(60:120))/mean(cres12(60:120))

print,'' 

print, mean(mres01(60:120))/mean(mres11(60:120))
print, mean(mres02(60:120))/mean(mres12(60:120))

print,'' 

print, mean(eres01(60:120))/mean(eres11(60:120))
print, mean(eres02(60:120))/mean(eres12(60:120))
end