par
# compute and plot Parallactic angle
# 
# Filippenko 1982 PASP 94 715
# by GA 25/02/2003 
# type "par" to start
               write standard  If you are at latitude=0, please move a bit...
               define emisfero ? {Northern emisphere "N", Southern emisphere "S"}
               

               if('$emisfero'=='N'){parnord} 
               if('$emisfero'=='S'){parsud}  
parnord

		

               define pi (acos(-1))
		
                define RA1 ? {RA hh}
                define RA2 ? {RA mm}
                define RA3 ? {RA ss}

                define RA ($RA1+($RA2+$RA3/60)/60) #RA decimal
                write standard RA=$RA            
 
                define ST1 ? {ST hh}
                define ST2 ? {ST mm}
                define ST3 ? {ST ss}

                define ST ($ST1+($ST2+$ST3/60)/60) #ST decimal
                write standard ST=$ST



               
                define h0 ($ST-$RA)
                if($h0<-12){define h0 ($h0+24)}
                if($h0>+12){define h0 ($h0-24)}

 
		define h0 ? {hour angle ? (decimal hours + west, - est)}
		#define long ? {observer longitudine ( + West, - Est; decimal degree) (-0.76 Ekar)}
                define fi ? {observer latitude (decimal degree, only + values!)(+45.83 Ekar, +42.26 Loiano, +33.36 Palomar, +25.76 TNG, +37.22 CAHA)}
 

                define delta1 ? {delta deg} #target declination (this value must be different from the latitude!!!)
                define delta2 ? {delta mm}
                define delta3 ? {delta ss}

                define delta ($delta1+($delta2+$delta3/60)/60) #delta decimal
                write standard delta=$delta


                
                define h (abs($h0))
		define hrad ($h*15/(180/$pi))
		define firad ($fi/(180/$pi))
                define deltarad ($delta/(180/$pi))
                define deltamin ? {Minimum declination value to be plotted? (-40)}
      		define deltamax ? {Maximum declination value to be plotted? (<90)}
     	        define step ? {Declination step? (10)}        

                
                erase
                limits 0 8 0 180
                ctype yellow
                box
                xlabel H.A.
                ylabel P.A.
                
                set xline=0,8,.1
                set yline1=xline*0
                set yline2=xline*0
                set yline3=xline*0
                define n 0
	        do hh=0,8,.1 {
		
                	
			define hrad ($hh*15/(180/$pi))
			define firad ($fi/(180/$pi))
              	        define deltarad ($delta/(180/$pi))
		
			define etarad (asin((sin($hrad)*cos($firad))/sqrt(1-(sin($firad)*sin($deltarad)+cos($firad)*cos($deltarad)*cos($hrad))**2)))

			define eta1 ($etarad*(180/$pi))
            	        define eta2 (180-$eta1)
		 	define eta3 (-180-$eta1)
			

                        set x=$hh
                        
                        set yline1[$n]=$eta1
                        set yline2[$n]=$eta2
                        set yline3[$n]=$eta3
                            
                       define n ($n+1)
	                }
                        
                ctype blue
                lweight 4
                if($delta<$fi){
                               set yline4= yline1
                               connect xline yline4
                               }
                
                if($delta>$fi){
                              set xlined=xline
                              set diff= yline2-yline1   
                              sort {diff xlined}
                              define  x1 (xlined[0])
                              set yline4= xline<$x1 ? yline2 : yline1 
                              connect xline  yline4
                               }
                if($delta==$fi){write standard *******************************************************
                               write standard CAUTION!!!!
                               write standard if Dec=$delta=Latitude=$fi
                               write standard ambigous angle
                               write standard *******************************************************
                               relocate (7000 20000)
                               ctype black
                               label Please enter declination different from the latitude!
                               }

              interp xline yline4 $h PARALLACTIC_ANGLE
        
              if($h0<0) {set PARALLACTIC_ANGLE=- PARALLACTIC_ANGLE}     

		write standard ************************************************
                write standard
		write standard latitude=$fi (decimal degree)
		write standard dec=$delta   (decimal degree)
                 write standard HOUR ANGLE=$h0
                               print {PARALLACTIC_ANGLE}
                 write standard WARNING: AFOSC slit is oriented  E-W (and not N-S as usual) so: 
                 #if(PARALLACTIC_ANGLE>180) {set AFOSC_PA=PARALLACTIC_ANGLE-180-90}  180<P.A<180
                 if(PARALLACTIC_ANGLE<180&&PARALLACTIC_ANGLE>0) {set AFOSC_PA=PARALLACTIC_ANGLE-90}  
                 if(PARALLACTIC_ANGLE<0&&PARALLACTIC_ANGLE>-180) {set AFOSC_PA=PARALLACTIC_ANGLE+180-90}  
                 #if(PARALLACTIC_ANGLE<-180){ set AFOSC_PA=PARALLACTIC_ANGLE-180-90}                     
                 print {AFOSC_PA}
		write standard ************************************************



                  set xlab=4
                  set delta22=$delta
                  set ynome=string(delta22)
                  ptype ynome
                  interp xline yline4 xlab  ylab
                  set ylab=ylab+5
                  points xlab ylab     


########################################### plot for different values of declination  ###################################

       lweight 1  

       do  delta2=$deltamin,$deltamax,$step { 
                if($delta2==$fi){define delta2 ($fi+0.00001)
                                write standard =======================================
                                write standard ATTENTION:
                                write standard if  declination = latitude
                                write standard not reliable values!!!!!!!!!!!!!!!!!!
                                write standard =======================================}

                define delta $delta2         
                set xline=0,8,.1
                set yline1=xline*0
                set yline2=xline*0
                set yline3=xline*0
                
                define n 0
	        do hh=0,8,.1 {
                	define hh (abs($hh))
			define hrad ($hh*15/(180/$pi))
			define firad ($fi/(180/$pi))
              	        define deltarad ($delta/(180/$pi))
		
			define etarad (asin((sin($hrad)*cos($firad))/sqrt(1-(sin($firad)*sin($deltarad)+cos($firad)*cos($deltarad)*cos($hrad))**2)))

			define eta1 ($etarad*(180/$pi))
            	        define eta2 (180-$eta1)
		 	define eta3 (-180-$eta1)
			

                        set x=$hh
                                                
                        set yline1[$n]=$eta1
                        set yline2[$n]=$eta2
                        set yline3[$n]=$eta3
                            
                        
                       define n ($n+1)
	                }


		

                ctype white
                if($delta<$fi){
                               set yline4= yline1
                               connect xline yline4
                               }
                
                if($delta>$fi){
                              set xlined=xline
                              set diff= yline2-yline1   
                              sort {diff xlined}
                              define  x1 (xlined[0])
                              set yline4= xline<$x1 ? yline2 : yline1 
                              connect xline  yline4
                               }
                if($delta==$fi){write standard *******************************************************
                               write standard CAUTION!!!!
                               write standard if Dec=$delta=Latitude=$fi
                               write standard ambigous angle
                               write standard *******************************************************
                               relocate (7000 20000)
                               ctype white
                               label Please enter declination different from the latitude!
                               }
                  
                  set xlab=2+(abs($delta2)/100)
                  set delta22=$delta2
                  set ynome=string(delta22)
                  ptype ynome
                  interp xline yline4 xlab  ylab
                  set ylab=ylab+5

                  points xlab ylab     

                         }       

                  write standard DONE

##################cambio emisfero

parsud

               define pi (acos(-1))
			
                define RA1 ? {RA hh}
                define RA2 ? {RA mm}
                define RA3 ? {RA ss}

                define RA ($RA1+($RA2+$RA3/60)/60) #RA decimal
                write standard RA=$RA            
 
                define ST1 ? {ST hh}
                define ST2 ? {ST mm}
                define ST3 ? {ST ss}

                define ST ($ST1+($ST2+$ST3/60)/60) #ST decimal
                write standard ST=$ST



               
                define h0 ($ST-$RA)
                if($h0<-12){define h0 ($h0+24)}
                if($h0>+12){define h0 ($h0-24)}

                
		
		#define long ? {observer longitudine ( + West, - Est; decimal degree)}
		define fi ? {observer latitude (decimal degree, only - values!) (-29.25 La Silla, -24.67 Paranal, -29.0 Las Campanas}
                

                define delta1 ? {delta deg} #target declination (this value must be different from the latitude!!!)
                define delta2 ? {delta mm}
                define delta3 ? {delta ss}

                define delta ($delta1+($delta2+$delta3/60)/60) #delta decimal
                write standard delta=$delta




                define h (abs($h0))
		define hrad ($h*15/(180/$pi))
		define firad ($fi/(180/$pi))
                define deltarad ($delta/(180/$pi))
                define deltamin ? {Minimum declination value to be plotted? (>-90)}
      		define deltamax ? {Maximum declination value to be plotted? (+40)}
     	        define step ? {Declination step? (10)}        

              
                
                erase
                limits 0 8 0 180
                ctype yellow
                box
                xlabel A.H.
                ylabel P.A.
                
                set xline=0,8,.1
                set yline1=xline*0
                set yline2=xline*0
                set yline3=xline*0
                define n 0
	        do hh=0,8,.1 {
		
                	
			define hrad ($hh*15/(180/$pi))
			define firad ($fi/(180/$pi))
              	        define deltarad ($delta/(180/$pi))
		
			define etarad (asin((sin($hrad)*cos($firad))/sqrt(1-(sin($firad)*sin($deltarad)+cos($firad)*cos($deltarad)*cos($hrad))**2)))

			define eta1 ($etarad*(180/$pi))
            	        define eta2 (180-$eta1)
		 	define eta3 (-180-$eta1)
			

                        set x=$hh
                        
                        set yline1[$n]=$eta1
                        set yline2[$n]=$eta2
                        set yline3[$n]=$eta3
                            
                       define n ($n+1)
	                }
                       #ptype 3 0 
                       #ctype red
                       #points xline yline1
                       #ctype green
                       #points xline yline2
                       #ctype cyan
                       #points xline yline3

                ctype blue
                lweight 4
                if($delta>$fi){ 
                               set yline4= yline2
                               connect xline yline4
                               }
                
                if($delta<$fi){
                              set xlined=xline
                              set diff= yline2-yline1   
                              sort {diff xlined}
                              define  x1 (xlined[0])
                              set yline4= xline<$x1 ? yline1 : yline2 
                              connect xline  yline4
                               }
                if($delta==$fi){write standard *******************************************************
                               write standard CAUTION!!!!
                               write standard if Dec=$delta=Latitude=$fi
                               write standard ambigous angle
                               write standard *******************************************************
                               relocate (7000 20000)
                               ctype white
                               label Please enter declination different from the latitude!
                               }

              interp xline yline4 $h PARALLACTIC_ANGLE
        
              if($h0<0) {set PARALLACTIC_ANGLE=- PARALLACTIC_ANGLE}     

		write standard ************************************************
                write standard
		write standard latitude=$fi (decimal degree)
		write standard dec=$delta   (decimal degree)
                write standard HOUR ANGLE=$h0
                               print {PARALLACTIC_ANGLE}
		write standard ************************************************



                  set xlab=4
                  set delta22=$delta
                  set ynome=string(delta22)
                  ptype ynome
                  interp xline yline4 xlab  ylab
                  set ylab=ylab+5
                  points xlab ylab     



########################################### general plot ###################################

       lweight 1  

       do  delta2=$deltamin,$deltamax,$step { 
                if($delta2==$fi){define delta2 ($fi+0.00001)
                                write standard =======================================
                                write standard ATTENTION:
                                write standard if  declination = latitude
                                write standard not reliable values!!!!!!!!!!!!!!!!!!
                                write standard =======================================}

                define delta $delta2         
                set xline=0,8,.1
                set yline1=xline*0
                set yline2=xline*0
                set yline3=xline*0
                
                define n 0
	        do hh=0,8,.1 {
                	define hh (abs($hh))
			define hrad ($hh*15/(180/$pi))
			define firad ($fi/(180/$pi))
              	        define deltarad ($delta/(180/$pi))
		
			define etarad (asin((sin($hrad)*cos($firad))/sqrt(1-(sin($firad)*sin($deltarad)+cos($firad)*cos($deltarad)*cos($hrad))**2)))

			define eta1 ($etarad*(180/$pi))
            	        define eta2 (180-$eta1)
		 	define eta3 (-180-$eta1)
			

                        set x=$hh
                                                
                        set yline1[$n]=$eta1
                        set yline2[$n]=$eta2
                        set yline3[$n]=$eta3
                            
                        
                       define n ($n+1)
	                }


		

                ctype white
                if($delta>$fi){
                               set yline4= yline2
                               connect xline yline4
                               }
                
                if($delta<$fi){
                              set xlined=xline
                              set diff= yline2-yline1   
                              sort {diff xlined}
                              define  x1 (xlined[0])
                              set yline4= xline<$x1 ? yline1 : yline2 
                              connect xline  yline4
                               }
                if($delta==$fi){write standard *******************************************************
                               write standard CAUTION!!!!
                               write standard if Dec=$delta=Latitude=$fi
                               write standard ambigous angle
                               write standard *******************************************************
                               relocate (7000 20000)
                               ctype white
                               label Please enter declination different from the latitude!
                               }
                  
                  set xlab=2+(abs($delta2)/100)
                  set delta22=$delta2
                  set ynome=string(delta22)
                  ptype ynome
                  interp xline yline4 xlab  ylab
                  set ylab=ylab+5

                  points xlab ylab     

                         }       
                  write standard DONE
