phase_compare_v2.R

phase_compare2<-function(GT,gal4){
  
  par(xpd=TRUE)
  par(mai=c(0.1,0.1,0.1,0.1))
  par(oma=c(2,0,1.5,0))
  
  gtcol=c("black","gray","black","gray")
  #gtcol=c(gray(0),gray(0.5),gray(.8),gray(.9))
  pty=c(16,16,8,8)
  #gtlty=c(1,1,3,3)
  gtlty=c(1,1,1,3)
  source('~/Dropbox/R/by_date2_V2.R')
  source('~/Dropbox/R/phaseplot.R')
  source('~/Dropbox/R/ccf_phase_plot2_V4.R')
  source('~/Dropbox/R/useful/ccf_phase_plot2_ebv2.R')
  library("circular")
  library("matrixStats")
  
  if(gal4=="tim"){cell=c(1:5)}
  if(gal4=="pdf"){cell=c(2,3)}
  if(gal4=="pdfrB"){cell=c(1,6,3)}
  if(gal4=="pdf80"){cell=c(4,5,1)}
  if(gal4=="pdfrB_pdf80"){cell=c(1,6)}
  if(gal4=="929"){cell=c(2)}
  lgd2=c("LNd","l-LNv","s-LNv","DN3","DN1","5th s-LNv")
  color=c("blue","red","orange","magenta","green","cyan")
  
  if(gal4=="eb"){
    cell=c(1,2)
    lgd2=c("evening","morning")
    color=c("blue","orange")
    ii=1
    }
  
  if(gal4=="TH"){
    cell=c(1,2,4,5)
    lgd2=c("PAL","PAM","PPM3","PPM1","PPL1")
    color=c("gold","orchid2","seagreen2","tomato","deepskyblue2")
    }
  
  if(gal4=="PPM3"){
    cell=c(1,2)
    lgd2=c("evening","morning")
    color=c("blue","orange")
    ii=5
  }
  if(is.numeric(gal4)){
    #cell=seq(1,gal4)
    cell=gal4
    cells=5}
  
  for (k in 1:length(cell)){
    if(gal4=="eb"||gal4=="PPM3"){ccfphase1=ccf_phase_plot2_eb(GT[1],ii=ii,plot="F")}
    else{ccfphase1=ccf_phase_plot2(GT[1],gal4,plot="F")}
    
    res=plot(ccfphase1[[k]],stack=F,bins=144,cex=1,col=gtcol[1],axes=F)
    #res=plot(pp1,stack=TRUE,bins=48,col="orange",pch=8,shrink=1.5,cex=1,axes = F)
    axis.circular(at=circular(seq(pi/4, 2*pi, pi/4))+pi/2, labels=seq(21,0,-3))
    
    arrows.circular(mean(ccfphase1[[k]],na.rm=TRUE),rho.circular(ccfphase1[[k]],na.rm=TRUE),
                    plot.info=res,lwd=6,col=gtcol[1])
    for(gt in 2:length(GT)){
      #ccfphase2=ccf_phase_plot2(GT[gt],gal4,plot="F")
      
      if(gal4=="eb"||gal4=="PPM3"){ccfphase2=ccf_phase_plot2_eb(GT[gt],ii=ii,plot="F")}
      else{ccfphase2=ccf_phase_plot2(GT[gt],gal4,plot="F")}
      
      points.circular(ccfphase2[[k]], col=gtcol[gt],plot.info=res,next.points=0.03*gt,cex=1,pch=pty[gt])
      
      
      arrows.circular(mean(ccfphase2[[k]],na.rm=TRUE),rho.circular(ccfphase2[[k]],na.rm=TRUE),
                      plot.info=res,lwd=6,col=gtcol[gt],lty=gtlty[gt])
      if(watson.williams.test(list(ccfphase1[[k]],ccfphase2[[k]]))$p.value<0.05){
        p1=mean(ccfphase1[[k]],na.rm=T)
        p2=mean(ccfphase2[[k]],na.rm=T)
        if(p1<0){p1=p1+24}
        if(p2<0){p2=p2+24}
        if(abs(p1-p2)>12){p2=p2-24}
        delay=1
        if(p1<p2){
          x=circular(seq(p2,p1,-0.01*delay),units="hours",template="clock24")
        }
        else{x=circular(seq(p1,p2,-0.01*delay),units="hours",template="clock24")}
         lines(x,rep((-.5-0.1*gt),length(x)),plot.info=res,lwd=3,col="black")
      }
      #mtext(text=sprintf("%s: p=%s",lgd2[cell[k]],signif(watson.williams.test(list(ccfphase1[[k]],ccfphase2[[k]]))$p.value,2)), side=1, line=0.3+gt-2,col=color[cell[k]] ,las=1, cex=3)
      #mtext(text=sprintf("%s: p=%s",nameGT[gt],signif(watson.williams.test(list(ccfphase1[[k]],ccfphase2[[k]]))$p.value,2)), side=1, line=gt-2,col=gtcol[gt] ,las=1, cex=2)
      
    }
    mtext(text=sprintf("%s: p=%s",lgd2[cell[k]],signif(watson.williams.test(list(ccfphase1[[k]],ccfphase2[[k]]))$p.value,2)), side=1, line=0.3,col=color[cell[k]] ,las=1, cex=3)
    #title(main=sprintf("%s",lgd2[cell[k]]),cex=1,col.main=color[k])
    #mtext(text=sprintf("%s",lgd2[cell[k]]), side=3, line=0,col=color[cell[k]] ,las=1, cex=3)
  }
}